function [pc5T, pc5W, pc5H] = precor55(Fun, t0, tlast, w0, hmax, hmin, TOL)
% %fc5T, pc5W, pc5H] = precor55(Fun, t0, tlast, w0, hmax, hmin, TOL)
%       Predictor - Corrector order 4 method based on Adams-Bashford
%	four step method and Adams-Moulton three step method with 
%       variable step size to solve  
%	    y' = Fun(t,y), on [t0, tlast] with y(t0)=w0 
%
% INPUT:
% Fun   - String containing name of function m-file (enclose in '  ')
% t0    - Initial value of t.
% tlast - Last value of t.
% w0    - Initial value.
% hmax  - maximum step size;  defaults to (tlast-t0)/10
% hmin  - minimum step size;  defaults to (tlast-t0)/200
% TOL   - local truncation tolerence; defaults to 1.e-6
%
% OUTPUT:
% pc5T  - time points (column-vector).
% pc5W  - Solution column-vector per pc5T-values.
% pc5H  - Column of step sizes h
% The result can be displayed by: plot(pc5T, pc5H).

if nargin < 4, 
 error('Not enough arguments, enter help precor55 for guidelines')
end;
if nargin < 5, 
 hmax=(tlast-t0)/10;
end;
if nargin < 7, 
 TOL = 1.e-6;
end;
if nargin < 6, 
 hmin=(tlast-t0)/200;
end;
% Initialization
pc5T = t0;
h = hmax;
pc5W = w0;
pc5H=0;
FLAG=1;
LAST=0;

[rkt,rkw]=rk4(Fun,t0,t0+3*h,w0,3);
NFLAG=1;
t=rkt(4)+h;

while FLAG==1
   k0 = feval(Fun,rkt(4),rkw(4));
   k1 = feval(Fun,rkt(3),rkw(3));
   k2 = feval(Fun,rkt(2),rkw(2));
   k3 = feval(Fun,rkt(1),rkw(1));
   wp=rkw(4)+(h/24)*(55*k0-59*k1+37*k2-9*k3);
   k=feval(Fun,t,wp);
   wc=rkw(4)+(h/24)*(9*k+19*k0-5*k1+k2);
   sigma=19*abs(wc-wp)/(270*h);

   if sigma <= TOL;
          i=length(pc5T);
       if NFLAG==1
          if (i==1 | LAST==1)
            pc5W=[pc5W;rkw(2:4);wc];pc5T=[pc5T;rkt(2:4);t];pc5H=[pc5H;h;h;h;h];
            rkt=[rkt(2:4);t];rkw=[rkw(2:4);wc];
          else
             pc5W=[pc5W(1:i-3);rkw(2:4);wc];pc5T=[pc5T(1:i-3);rkt(2:4);t];
             pc5H=[pc5H(1:i-3);h;h;h;h];
             rkt=[rkt(2:4);t];rkw=[rkw(2:4);wc];
          end;
       else
          pc5W=[pc5W;wc];pc5T=[pc5T;t];pc5H=[pc5H;h];
          rkt=[pc5T(i-2:i);t];rkw=[pc5W(i-2:i);wc];
       end;
       if LAST==1 
                 FLAG=0;
       else
             NFLAG=0;
          if  (sigma <=.1*TOL | t +h > tlast)
                   q=(TOL/(2*sigma))^(1/4);
                      if q > 4
                          h=4*h;
                      else
                          h=q*h;
                      end;
                      if h > hmax
                         h=hmax;
                      end;
                      if t+4*h > tlast
                         h=(tlast - t)/4;
                         LAST=1;
                         [rkt,rkw]=rk4(Fun,t,t+3*h,rkw(4),3);
                         t=rkt(4); NFLAG=1;
                      end;
          end;
       end;  
   else
        q=(TOL/(2*sigma))^(1/4);LAST=0;
        if q < 0.1
                h=0.1*h
        else
                h=q*h;
        end;
        if h < hmin
                FLAG=0;
                error('hmin violated');
        else
                [rkt,rkw]=rk4(Fun,rkt(1),rkt(1)+3*h,rkw(1),3);
                t=rkt(4);NFLAG=1;
        end;
   end;
  t=t+h;
end;
