function [t,x]=pred_corr(field,t0,T,x0,h)
%PRED_CORR Implements the PECECE Predictor/Corrector algorithm using
%          Adams-Bashfort 2 (Predictor) and Adams-Moulton 2 (Corrector)
%
% Parameters:
%  field  -- Right hand side function of ODE system: x'=f(t,x)
%  t0     -- Initial time
%  T      -- End time (T > t0)
%  x0     -- Initial value
%  h      -- Size of time step (h <= T-t0)
%
% Outputs:
%  t  -- [t0; t-0+h, t0+2*h; ...; t0+i*h; ...]
%  x  -- Column vector containing numerical solution at each time step

q = 1;
n = ceil((T-t0)/h);

t = t0+h*(0:n).';
x = ones(n+1,length(x0));
[t(1:(q+1)), x(1:(q+1),:)] = rk_classical(field, t0, t0+q*h, x0, h);

f1 = feval(field, t(q), x(q,:).');
f0 = feval(field, t(q+1), x(q+1,:).');
for i=(q+1):n
    x0 = x(i,:).'+h*(3*f0-f1)/2;            % P - ab2
    f_new = feval(field, t(i,:)+h, x0);     % E
    x0 = x(i,:).' + h*(5*f_new+8*f0-f1)/12; % C - am2
    f_new = feval(field, t(i,:)+h, x0);     % E
    x0 = x(i,:).' + h*(5*f_new+8*f0-f1)/12; % C - am2
    f_new = feval(field, t(i,:)+h, x0);     % E
    
    x(i+1,:) = x0;
    f1 = f0;
    f0=f_new;
end
