% M-File grxxm103.m
% Steifes, nichtlineares AWP
% Runge-Kutta-Verfahren
% Steifigkeitsgrad über Jacobi-Matrix J(t)

clear;

f1=inline('-0.1*y1+100*y2*y3','y1','y2','y3');
f2=inline('0.1*y1-100*y2*y3-500*y2^2','y1','y2','y3');
f3=inline('500*y2^2-0.5*y3','y2','y3');

% Eigenwerte von J(0)
J=[-0.1 50 200;0.1 -2050 -200;0 2000 -0.5];
lam=eig(J);
slam=sort(lam);

y(1,1:8)=[0,4,2,0.5,slam(1),slam(2),slam(3),slam(1)/slam(3)];

% 25 Runge-Kutta-Schritte mit h=0.0002
% 995 Runge-Kutta-Schritte mit h=0.001
h=0.0002;


for k=1:1020
  if k<26 
      h=0.0002;
  else
      h=0.001;
  end;    
    y(k+1,1)=h*k;
    k1=[f1(y(k,2),y(k,3),y(k,4)),...
        f2(y(k,2),y(k,3),y(k,4)),...
        f3(y(k,3),y(k,4))];
    k2=[f1(y(k,2)+h*k1(1)/2,y(k,3)+h*k1(2)/2,y(k,4)+h*k1(3)/2),...
        f2(y(k,2)+h*k1(1)/2,y(k,3)+h*k1(2)/2,y(k,4)+h*k1(3)/2),...
        f3(y(k,3)+h*k1(2)/2,y(k,4)+h*k1(3)/2)];
    k3=[f1(y(k,2)+h*k2(1)/2,y(k,3)+h*k2(2)/2,y(k,4)+h*k2(3)/2),...
        f2(y(k,2)+h*k2(1)/2,y(k,3)+h*k2(2)/2,y(k,4)+h*k2(3)/2),...
        f3(y(k,3)+h*k2(2)/2,y(k,4)+h*k2(3)/2)];
    k4=[f1(y(k,2)+h*k3(1),y(k,3)+h*k3(2),y(k,4)+h*k3(3)),...
        f2(y(k,2)+h*k3(1),y(k,3)+h*k3(2),y(k,4)+h*k3(3)),...
        f3(y(k,3)+h*k3(2),y(k,4)+h*k3(3))];
    w=(k1+2*k2+2*k3+k4)*h/6;
        y(k+1,2)=y(k,2)+w(1);
        y(k+1,3)=y(k,3)+w(2);
        y(k+1,4)=y(k,4)+w(3);
        J=[-0.1 100*y(k+1,4) 100*y(k+1,3);...
            0.1 -100*y(k+1,4)-1000*y(k+1,3) -100*y(k+1,3);
            0 1000*y(k+1,3) -0.5];
        lam=eig(J);
        slam=sort(lam);
        y(k+1,5)=slam(1);
        y(k+1,6)=slam(2);
        y(k+1,7)=slam(3);
        y(k+1,8)=slam(1)/slam(3);
end;

% Tabelle 
disp('------------------------------------------------------------------------------------------------');
disp('       Runge - Kutta - Verfahren für das steife, nichtlineares AWP aus Aufgabe 3');
disp('------------------------------------------------------------------------------------------------');
disp(' t_n      y_1n     y_2n     y_3n      lam_1n        lam_2n        lam_3n         S_n    ');
disp('------------------------------------------------------------------------------------------------');

for s=1:5:50
    fprintf('%1.4f %8.4f %8.4f %8.4f %10.6e %10.6e %13.6e %13.6e',...
        y(s,1),y(s,2),y(s,3),y(s,4),y(s,5),y(s,6),y(s,7),y(s,8));
    fprintf('\n');
end;
for s=51:50:1001
    fprintf('%1.4f %8.4f %8.4f %8.4f %10.6e %10.6e %13.6e %13.6e',...
        y(s,1),y(s,2),y(s,3),y(s,4),y(s,5),y(s,6),y(s,7),y(s,8));
    fprintf('\n');
end;
