
function [x,t]=trapzPC(f,t0,T,x0,N=20,p=1)
%predictor corrector based on forward Heun -trapezoid  schemes
%base on p step of Bacnach iteration for trapezoid i.e.
%for given x(n) x^0=x(n)+0.5*h*(f(n) +f(x(n)+h*f(n),t(n+1)) - one step of forward Euler
% x(n+1)=x^p for x^(l+1)=x^l+0.5*h*(f(n)+f(x^l,t(n+1)))
%EXAMPLE
% [x,t]=trapzPC(f=@(x,t) -x,t0=0,T=1,x0=1,N=20,p=1)
  h=(T-t0)/N;
  m=length(x0);
  t=linspace(t0,T,N+1);
  x=zeros(N+1,m)';
  x(:,1)=x0;
  for k=1:N,
     fx=f(x(:,k),t(k));
     xx=x(:,k)+h*fx; %forward Euler
     xx=x(:,k)+0.5*h*(fx+ f(xx,t(k))); %Heun 
     for l=1:p,
      %x(n+1)=x(n)+h*f(x(n+1),t(n+1))
      xx=x(:,k)+0.5*h*(fx+f(xx,t(k))); %Banach iteration for backward Euler
     end
     x(:,k+1)=xx;
  end   
  x=x'; %for conformity to the octave black box ode solver
end