function [pc4T, pc4W] = precor54(Fun, t0, tlast, w0, N)
%%function [pc4T, pc4W] = precor54(Fun, t0, tlast, w0, N)
%       Predictor - Corrector order 4 method based on Adams-Bashford
%	four step method and Adams-Moulton three step method to solve
%	    y' = Fun(t,y), on [t0, tlast] with y(t0)=w0 with N+1 points
% INPUT:
% Fun   - String containing name of function m-file (enclose in '  ')
% t0    - Initial value of t.
% tlast - Last value of t.
% w0    - Initial value.
% N     - The number of steps
%
% OUTPUT:
% pc4T  - time points (column-vector).
% pc4W  - Solution column-vector per pc4T-values.
%
% The result can be displayed by: plot(pc4T, pc4W).

if nargin < 5, 
 error('Not enough arguments, enter help precor54 for help')
end;
% Initialization
pc4T = t0;
h = (tlast - t0)/N;
pc4W = w0;
%% RK of order 4 to find first few steps
for i=1:3;
 k1= h*feval(Fun,pc4T(i),pc4W(i));  
 k2= h*feval(Fun,pc4T(i)+(h/2),pc4W(i)+(k1/2));  
 k3= h*feval(Fun,pc4T(i)+(h/2),pc4W(i)+(k2/2));  
 k4= h*feval(Fun,pc4T(i)+h,pc4W(i)+k3);  
 w = pc4W(i) + (k1+2*k2+2*k3+k4)/6; 
 t = pc4T(i) + h; clear k1 k2 k3 k4;
 pc4T=[pc4T;t]; pc4W = [pc4W;w];
end;

% The main loop
   for i=4:N;
      k_i = feval(Fun,pc4T(i),pc4W(i));
      k_i1 = feval(Fun,pc4T(i-1),pc4W(i-1));
      k_i2 = feval(Fun,pc4T(i-2),pc4W(i-2));
      k_i3 = feval(Fun,pc4T(i-3),pc4W(i-3));
      w=pc4W(i)+(h/24)*(55*k_i-59*k_i1+37*k_i2-9*k_i3);
         k=feval(Fun,pc4T(i)+h,w);
      w=pc4W(i)+(h/24)*(9*k+19*k_i-5*k_i1+k_i2);
      t=t+h;
      pc4T=[pc4T;t]; pc4W = [pc4W;w];
   end;

