function z = bdf ( k, maxphi, points );
% BDF
%
% 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 von BDF
%
% Parameter
%  k      : Grad der BDF (k<=7)
%  maxphi : Intervallende bez. phi in Vielfachen von 2*pi (Anzahl der Umlaeufe um den Koordinatenursprung)
%  points : Anzahl der Berechnungspunkte

% Verfahrenskoeffizienten
d=zeros(8,8);
d(:,1)=cos(zeros(8,1));
for l1=1:7,
  for l2=2:l1+1,
    d(l1+1,l2)=d(l1,l2)-d(l1,l2-1);
  end;
end;
coeff=zeros(1,8);
for l=1:7,
  coeff(l+1)=1/l;
  bdfk(l,:)=coeff*d;
end;

% Initialisierung   
phi=linspace(0, 2*maxphi*pi, points);

% Kurvenverfolgung
for point=1:points,     
  xi=exp(i*phi(point));
  r=bdfk(k, k+1);
  for l=1:k,
    r=r/xi+bdfk(k, k-l+1);
  end;  
  z(point)=r;                            
end;     

