% Oblast stability AB2 metody
% y_{i+1}=y_i+h(3/2y'_i-1/2y'_{i-1})
% kořeny: r0=1/2+3/4z+sqrt(1+z+9/4z^2), r1=1/2+3/4z-sqrt(1+z+9/4z^2),

r0=@(z)abs(1/2+3/4*z+0.5*sqrt(1+z+9/4*z.^2));
r1=@(z)abs(1/2+3/4*z-0.5*sqrt(1+z+9/4*z.^2));
x=linspace(-2,2,41);
y=linspace(-2,2,41);
[X,Y]=meshgrid(x,y);
Z=X+i*Y;
R0=r0(Z);
mesh(x,y,R0)
contour(x,y,R0,[1 1])

%
%

R1=r1(Z);
mesh(x,y,R1)
contour(x,y,R1,[1 1])

%
%

P=R1./R0;

mesh(x,y,P)
contour(x,y,P,[1 1])

%
% BDF2 metoda

r0=@(z)abs(0.5*(4/3+2/3*sqrt(1+2*z))./(1-2/3*z));
r1=@(z)abs(0.5*(4/3-2/3*sqrt(1+2*z))./(1-2/3*z));
x=linspace(-3,1,41);
y=linspace(-2,2,41);
[X,Y]=meshgrid(x,y);
Z=X+i*Y;
R0=r0(Z);
mesh(x,y,R0)
contour(x,y,R0,[1 1])

%
%

R1=r1(Z);
mesh(x,y,R1)
[C,h]=contour(x,y,R1);
clabel(C,h)

%

x=linspace(-3,3,61);
y=linspace(-2,2,41);
[X,Y]=meshgrid(x,y);
Z=X+i*Y;
R0=r0(Z);
R1=r1(Z);
P=R1./R0;

mesh(x,y,P)
[C,h]=contour(x,y,P);
clabel(C,h)

%
%

R0=r0(Z);
mesh(x,y,R0)
[C,h]=contour(x,y,R0,[1 2 3 4 5]);
clabel(C,h)