% grxxm083.m
% Räuber - Beute - Modell nach V.Volterra
% implizites Zweischrittverfahren (Milne-Simpson-Verfahren)
% u_{n+1} = u_{n-1} + h/3 ( f_{n+1} + 4 f_n + f_{n-1} )

clear;

%% Grafikfenster (nur MATLAB / OCTAVE: ... Kommentar)
%% ... eine huebsche GUI-Oberflaeche fuer MATLAB ...
% screen=get(0,'ScreenSize');
% screenplot=[50 50 screen(3)*0.8 screen(4)*0.8];
% figure('Position',screenplot,'color',[1 1 0.77]);
%% Schliessen Grafikfenster (nur MATLAB / OCTAVE: ... Kommentar)
% 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');

%% ab hier: beide Systeme arbeiten
     
f1=inline('-a*x+b*x*y','x','y','a','b');
f2=inline('c*y-d*x*y','x','y','c','d');
     
% Modelldaten
a=0.2;
b=0.03;
c=0.1;
d=0.03;

% Anfangswerte
x0=6;
y0=15;

%x0=c/d;
%y0=a/b;

%  Abbruchschranke für die Iteration
delta=1.0e-8;  

% Schrittweite h
T=70;
n=7000;
h=T/n;

r(1,1:5)=[0 x0 y0 f1(x0,y0,a,b) f2(x0,y0,c,d)];

% Anlaufrechnung mit R-K-Verfahren
k1=[f1(x0,y0,a,b);f2(x0,y0,c,d)];
k2=[f1(x0+k1(1)*h/2,y0+k1(2)*h/2,a,b);f2(x0+k1(1)*h/2,y0+k1(2)*h/2,c,d)];
k3=[f1(x0+k2(1)*h/2,y0+k2(2)*h/2,a,b);f2(x0+k2(1)*h/2,y0+k2(2)*h/2,c,d)];
k4=[f1(x0+k3(1)*h,y0+k3(2)*h,a,b);f2(x0+k3(1)*h,y0+k3(2)*h,c,d)];
dx=[x0;y0]+(k1+2*k2+2*k3+k4)*h/6;
r(2,1:5)=[h dx(1) dx(2) f1(dx(1),dx(2),a,b) f2(dx(1),dx(2),c,d)];

% implizites Zweischrittverfahren
for k=2:n
    r(k+1,1)=k*h;
      u=[r(k-1,2);r(k-1,3)]+([r(k-1,4);r(k-1,5)]+4*[r(k,4);r(k,5)])*h/3;
    v=[r(k,2);r(k,3)];
    s=1;
  while s>delta
    w=u+[-a*v(1)+b*v(1)*v(2);c*v(2)-d*v(1)*v(2)]*h/3;
    s=norm(w-v,2);
    v=w;
  end;
    r(k+1,2)=w(1);
    r(k+1,3)=w(2);
    z=[-a*w(1)+b*w(1)*w(2);c*w(2)-d*w(1)*w(2)];
    r(k+1,4)=z(1);
    r(k+1,5)=z(2);
end;  

% Plots
subplot(2,2,1);
plot(r(:,1),r(:,3),'b');
xlabel('Zeit  t','FontSize',12);
ylabel('y(t)','FontSize',12);
title('Beutespezies  y(t)','FontSize',15);

subplot(2,2,3);
plot(r(:,1),r(:,2),'r');
xlabel('Zeit  t','FontSize',12);
ylabel('x(t)','FontSize',12);
title('Raeuberspezies  x(t)','FontSize',15);


subplot(2,2,2);
plot(r(:,1),r(:,2),'r',r(:,1),r(:,3),'b');
xlabel('Zeit  t','FontSize',12);
ylabel('x(t),  y(t)','FontSize',12);
title('Raeuber-  und  Beutespezies','FontSize',15);
legend('Raeuber','Beute');

subplot(2,2,4);
plot(r(:,2),r(:,3),'m');
hold on;
plot(x0,y0,'ok');
xlabel('Raeuberspezies  x(t)','FontSize',12);
ylabel('Beutespezies y(t)','FontSize',12);
title('x-y-Diagramm','FontSize',15);
