function  z = explrk ( method, maxphi, points, z0 );
% EXPLRK
%
% This m-file is part of the material of the lecture
%   Scientific Computing and Software Tools
% at the Math Department of Halle University, Germany.
% Summer term 2004
%
% Author :      Prof. Dr. M. Arnold, arnold@mathematik.uni-halle.de
% Version of :  May 11, 2004
%
% Berechnung und graphische Darstellung des Stabilitaetsbereichs expliziter
% Runge-Kutta-Verfahren
%
% Parameter
%  method : explizites Runge-Kutta-Verfahren
%           1 - Euler
%           2 - DoPri 5(4)
%  maxphi : Intervallende bez. phi in Vielfachen von 2*pi (Anzahl der Umlaeufe um den Koordinatenursprung)
%  points : Anzahl der Berechnungspunkte      
%  z0     : Anfangspunkt fuer die Berechnung (geschaetzt) , z.B. z0=0

% Verfahrenskoeffizienten
aeuler=[0];
beuler=[1];
adopri= ...
  [      0           0          0         0         0        0   ,
        1/5          0          0         0         0        0   ,
        3/40        9/40        0         0         0        0   ,
       44/45      -56/15      32/9        0         0        0   ,
    19372/6561 -25360/2187 64448/6561 -212/729      0        0   ,
     9017/3168   -355/33   46732/5247   49/176 -5103/18656   0   ];
bdopri= ...
  [    35/384        0       500/1113  125/192 -2187/6784  11/84 ];
if method==1
  a=aeuler;
  b=beuler;
elseif method==2
  a=adopri;
  b=bdopri;
end;

% Polynomkoeffizienten (Stabilitaetsfunktion)
s=length(b);
e=cos(zeros(s, 1));
for k=1:s,
  coeff(k)=b*e;
  b=b*a;
end;     

% Initialisierung   
phi=linspace(0, 2*maxphi*pi, points);
zneu=z0;
eps=1e-10;

% Korrektor
diff=2*eps;
while diff>eps,
  r=coeff(s);
  for k=1:s-1,
    r=r*zneu+coeff(s-k);
  end;  
  r=r*zneu+1;
  rprime=s*coeff(s);
  for k=1:s-1,
    rprime=rprime*zneu+(s-k)*coeff(s-k);
  end;
  zaux=zneu-(r-1)/rprime;
  diff=abs(zaux-zneu);
  zneu=zaux;
end;  
z(1)=zneu;                            
% Kurvenverfolgung
for point=2:points,
% Praediktor
  zneu=z(point-1);
  rprime=s*coeff(s);
  for k=1:s-1,
    rprime=rprime*zneu+(s-k)*coeff(s-k);
  end;
  zneu=zneu+(phi(point)-phi(point-1))*i*exp(i*phi(point-1))/rprime; 
% Korrektor
  diff=2*eps;
  while diff>eps,
    r=coeff(s);
    for k=1:s-1,
      r=r*zneu+coeff(s-k);
    end;  
    r=r*zneu+1;
    rprime=s*coeff(s);
    for k=1:s-1,
      rprime=rprime*zneu+(s-k)*coeff(s-k);
    end;
    zaux=zneu-(r-exp(i*phi(point)))/rprime;
    diff=abs(zaux-zneu);
    zneu=zaux;
  end;  
  z(point)=zneu;                            
end;     
