clear all;
close all;
clc;

epsilon = [10.^[-8:-2],0.05:0.05:5];
Nrho = 120;
Nphi = 3600;

% w = z^2 - 3; proof of lim z->2 (w) = 1 

delta = zeros(length(epsilon),1);

figure;

for kk = 1 : length(epsilon)

    rho = epsilon(kk)/(2*Nrho): epsilon(kk)/Nrho : epsilon(kk) - epsilon(kk)/(2*Nrho);
    phi = pi/Nphi : 2*pi/Nphi : 2*pi - pi/Nphi;
    
    u = zeros(Nphi*Nrho,1); % Real part of w 
    v = zeros(Nphi*Nrho,1); % Imag part of w 
    pp = 1;
    for mm = 1 : Nphi
        for nn = 1 : Nrho
            u(pp) = 1 + rho(nn)*cos(phi(mm));
            v(pp) = rho(nn)*sin(phi(mm));
            pp = pp + 1;
        end
    end
    w = u + 1i*v;
    
    z = sqrt(w+3);
    
    x = real(z); y = imag(z);
    
    delta(kk) = max(abs(z-2));
    
    subplot(1,2,1);
    plot(1 + epsilon(kk)*cos([phi,pi/Nphi]), epsilon(kk)*sin([phi,pi/Nphi]), 'r-', u,v,'r*'); 
    title('\omega = z^2 - 3 = u + iv plane','FontSize',15);
    xlabel('u','FontSize',15);
    ylabel('v','FontSize',15);
    xlim([-7,7]);
    ylim([-7,7]);
%     annotation('arrow',[0.5,0.5],[0.2,0.8],'String','\epsilon ');
%     drawnow;
    subplot(1,2,2);
    plot(2+delta(kk)*cos([phi,pi/Nphi]),delta(kk)*sin([phi,pi/Nphi]),'k-',x,y,'k*');
    title('z = x + iy plane','FontSize',15);
    xlabel('x','FontSize',15);
    ylabel('y','FontSize',15);    
%     annotation('arrow',[2,0],[2,delta(kk)],'String',sprintf('\delta = %g ',delta(kk)));
    xlim([-7,7]);
    ylim([-7,7]);
    drawnow;
    hold off;
end

%delta = delta(epsilon) graph
figure;
plot(epsilon,delta,'b-');
title('There is a \delta(\epsilon)>0 for each \epsilon>0 !!!','FontSize',15);
xlabel('\epsilon','FontSize',15);
ylabel('\delta = \delta(\epsilon)','FontSize',15);
