function y = cauchy_integral(x)

L = length(x);
y = zeros(1,L);

N = 3000;
R = max(x(1),x(L))+1;
theta  = linspace(0,2*pi,N);
dtheta = theta(2)-theta(1); 
zeta   = R*exp(i*theta);
dzeta  = R*i*exp(i*theta)*dtheta;

% Cauchy integral formula
%------------------------------------------------------------
% Plot the (D-1)-th derivative of f(z)
D = 2;

for j = 1:L
    cauchy = 0;
    for k = 1:N
        term   = f(zeta(k))/(zeta(k)-x(j))^D*dzeta(k);
        cauchy = cauchy + term; 
    end
    y(j) = 1/(2*pi*i)*factorial(D-1)*cauchy;
end

% Results
%------------------------------------------------------------
plot(x,real(y),'r');
hold on;
y2 = zeros(1,L);
for j = 1:L
   y2(j) = f(x(j));
end
plot(x,y2,'b');
hold off;
xlabel('x');
ylabel('y');
title('Cauchy representation of f(z) over the real line (red) vs. original (blue)');
grid on;

% Function
%------------------------------------------------------------
function w = f(z)
w = z^3+3*z^2-4;