% Runge-Kutta

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

%
% Ralstonova 3.řádu

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

%
% TRX2

h=0.1;
k1=@(lambda,a11,a21,a22,a31,a32,a33) lambda./(1-lambda.*h.*a11);
k2=@(lambda,a11,a21,a22,a31,a32,a33) (lambda+lambda.*h.*a21.*k1(lambda,a11,a21,a22,a31,a32,a33))./(1-lambda.*h.*a22);
k3=@(lambda,a11,a21,a22,a31,a32,a33) (lambda+lambda.*h.*a31.*k1(lambda,a11,a21,a22,a31,a32,a33) ...
    +lambda.*h.*a32.*k2(lambda,a11,a21,a22,a31,a32,a33))./(1-lambda.*h.*a33);
f=@(lambda,a11,a21,a22,a31,a32,a33,b1,b2,b3) abs(1+h.*(b1.*k1(lambda,a11,a21,a22,a31,a32,a33)+ ...
    b2.*k2(lambda,a11,a21,a22,a31,a32,a33)+b3.*k3(lambda,a11,a21,a22,a31,a32,a33)));

%

a=-30:0.1:10;
b=-40:0.1:40;
a11=0;a21=1/4;a22=1/4;a31=1/4;a32=1/2;a33=1/4;b1=1/4;b2=1/2;b3=1/4;
my_f=@(lambda) f(lambda,a11,a21,a22,a31,a32,a33,b1,b2,b3)
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