% M-File grxxm093.m

% Effekte der Null-Instabilität
% AWP:   y' = (5-y)*y,  y(0)=1,
% Lösung  y(t) = 1/(1+4*exp(-5*t))

clear;

% Grafikfenster (nur MATLAB)
screen=get(0,'ScreenSize');
screenplot=[50 50 screen(3)*0.8 screen(4)*0.8];
figure('Position',screenplot,...
          'color',[1 1 0.77]);

clf;

% Schliessen Grafikfenster (nur MATLAB)
schliessx=(screenplot(3)-screenplot(1))/2;
schliess=[schliessx-60 10 120 30];
uicontrol(gcf,...
         'Style','Push',... 
         'Position',schliess,...
         'String','Schliessen',...
         'Backgroundcolor',[0.5 0.93 0.9],...
         'callback','close');
     
f=inline('(5-y)*y','y');
l=inline('5/(1+4*exp(-5*t))','t');

% Scrittweite h
T=2;
%n=20; %%%%   Auswahl von  "n"  %%%%
n=40;
%n=80;
h=T/n;

% Parameter "a" zur Beeinflussung der Nullstabilität
% -1 <= a < 1  Stabilitaetsbetreich
%%%%%%    Auswahl von a    %%%%%%
a=-1.1;
%a=0.4;
%a=1;
%a=1.1;
%a=0.4;
%a=-0.9;
%a=-1.0;
%a=-1.1;

y(1,1:5)=[0 1 4 1 0];

% Anlaufrechnung mit R-K-Verfahren
y(2,1)=h;
   k1=f(y(1,2));
   k2=f(y(1,2)+k1*h/2);
   k3=f(y(1,2)+k2*h/2);
   k4=f(y(1,2)+k3*h);
y(2,2)=y(1,2)+(k1+2*k2+2*k3+k4)*h/6;
y(2,3)=f(y(2,2));
y(2,4)=l(y(2,1));
y(2,5)=y(2,2)-y(2,4);

% Zweischrittformel (*)
for k=2:n
    y(k+1,1)=k*h;
      z=((3-a)*y(k,3)-(1+a)*y(k-1,3))*h/2;
    y(k+1,2)=(1+a)*y(k,2)-a*y(k-1,2)+z;
    y(k+1,3)=f(y(k+1,2));
    y(k+1,4)=l(y(k+1,1));
    y(k+1,5)=y(k+1,2)-y(k+1,4);
end;

% Plot: Loesung und Naeherung (abhaengig von a)
plot(y(:,1),y(:,2),'r',y(:,1),y(:,4),'b');
hnum=num2str(h);
anum=num2str(a);
str1=['Zweischrittverfahren:   Parameter  a = ' anum ',' ];
str2=['   Schrittweite  h = ' hnum];
title([str1 str2],'FontSize',15);
xlabel('t - Achse');
ylabel('y - Achse');
legend('Zweischrittverfahren','exakte Loesung',4);
