function [t,y,errest] = predcorr(f,tint,y0,h)
% [t,y,errest] = predcorr(f,tint,y0,h)
% Simple implementation of a predictor-corrector code 
% 
%   f : function handle to the right hand side of y'=f(t,y)
%   tint : Integration interval, tint = [t0,tend]
%   y0 : Initial values, y0 = y(t0) 
%   h  : Stepsize



% Some constants:
NITER = 2;                   % Number of iterations, your choice.

% Some initialization
t0 = tint(1);
tend = tint(2);
NSTEPS = floor((tend-t0)/h); % Number of steps needed.

t = t0;
y = y0;
F = feval(f,t0,y0);
errest = 0;

% Computation of y1, by improved Euler
k1 = feval(f,t0,y0);
k2 = feval(f,t0+h,y0+h*k1);
ytmp = y0 + 0.5*h*(k1+k2);
ftmp = feval(f,t0+h,ytmp);
erresttmp = 0;

y = [y ytmp];  % Store the solutions
t = [t t0+h];
F = [F ftmp];
errest = [errest erresttmp];

% Main loop:
for m=2:NSTEPS
  ytmp = y(:,m)+h*(1.5*F(:,m)-0.5*F(:,m-1));   % AB2
  ypred = ytmp;
  ttmp = t(m)+h;
  for i=1:NITER                 % PC-iterations
      ytmp = y(:,m) + 0.5*h*(F(:,m)+feval(f,ttmp,ytmp)); %(AB1)
  end
  ftmp = feval(f,ttmp,ytmp);
  erresttmp = -(ytmp-ypred)/6;  % Error estimation, by Milnes device. 

  y = [y ytmp];  % Store the solutions
  t = [t ttmp];
  F = [F ftmp];
  errest = [errest norm(erresttmp,2)];  
end
% End of main loop. 









