% Zona de estabilidad metodos numéricos RK3 stab_RK3.m, EQM
% metodo de solucion: Newton-Raphson
% datos: n=numero de iteraciones
%        x=valor inicial de la abscisa
%        y=valor inicial de la ordenada
%        h=paso de integracion
% nota: x=xh, y=yh, rev. 18enr'25
%----------------------------------------------------------------------
clear
close all
format short e
load xx4 % RK4
load xx3 % RK3
x=0.0001; yi=0; yf=-1.73; dy=-.01; xx=[]; % yf=-2.93, -2.37
for j=1:2
    if (j==2) yi=-1.73; yf=0; dy=.01; x=-1.001; end % valor inicial
    for i=yi:dy:yf
        y=i;
        loop=0;
        while (1)
            loop=loop+1;
            if (loop>100)
                error('no converge con 100 iteraciones');
            end
            Re=1+x+(1/2)*(x^2-y^2); %+(1/6)*(x^3-3*x*y^2);
                %(1/24)*(x^4-6*x^2*y^2+y^4);
            Im=y+(1/2)*(2*x*y); %+(1/6)*(-y^3+3*x^2*y);
                %(1/24)*(4*x^3*y-4*x*y^3);
            f=Re^2+Im^2-1;
            dRe=1+x; %+(1/2)*x^2-(1/2)*y^2; %+(1/6)*x^3-(1/2)*x*y^2;
            dIm=y; % +x*y; %+(1/2)*x^2*y-(1/6)*y^3;
            df=2*Re*dRe+2*Im*dIm;
            if(abs(f)<1e-10) break; end
            x=x-f/df;
        end
        xx=[xx;x y f loop];
    end
end
plot(xx4(:,1),xx4(:,2),xx4(:,1),-xx4(:,2)),grid,axis('equal',[-4 1 -3 3])
hold on
plot(xx3(:,1),xx3(:,2),xx3(:,1),-xx3(:,2))
plot(xx(:,1),xx(:,2),xx(:,1),-xx(:,2))
%title('Zona estabilidad metodo Runge-Kutta 4-3-2-1')
xlabel('Re(\lambdah)'); ylabel('Im(\lambdah)');
x1=0:-.01:-2;
y1=sqrt(1-(1+x1).^2);
plot(x1,y1,x1,-y1)
s=strvcat('4','3','2','1');
xx=[.3,-.2,-.5,-.9]'; yy=[2.6,2.3,1.7,1.2]';
text(xx,yy,s)
hold off