% Explicitní Eulerova

a=-30:0.1:30;
b=-30:0.1:30;
h=0.1;
f=@(lambda) abs(1+h.*lambda);
stability(a,b,f)

%
% Implicitní Eulerova

a=-30:0.1:30;
b=-30:0.1:30;
h=0.1;
f=@(lambda) abs(1-h.*lambda);
stability(a,b,f)

%
% Lichoběžníková

a=-40:0.1:10;
b=-30:0.1:30;
h=0.1;
f=@(lambda) abs((1+1/2.*h.*lambda)./(1-1/2.*h.*lambda));
stability(a,b,f)

%
% Runge-Kutta

a=-40:0.1:10;
b=-30:0.1:30;
h=0.1;
f=@(lambda,b1,b2,a21) abs(1+h*b1*lambda+h*b2*lambda+h.^2*b2*a21*lambda.^2);

%
% Modifikovaná EM

b1=0;b2=1;a21=1/2;
my_f=@(lambda) f(lambda,b1,b2,a21)
stability(a,b,my_f)

%
% Ralsnova

b1=1/4;b2=3/4;a21=2/3;
my_f=@(lambda) f(lambda,b1,b2,a21)
stability(a,b,my_f)

%
% Heunova

b1=1/2;b2=1/2;a21=1;
my_f=@(lambda) f(lambda,b1,b2,a21)
stability(a,b,my_f)

%

function []=stability(a,b,f)
[A,B]=meshgrid(a,b);
all_f=f(A+B*1i);
contour(a,b,all_f,10)
hold on
contour(a,b,all_f,[1 1],'r')
hold off
end