function [t,y,trunc,mil]=abmpc(a,b,eta,N,iex,k,m,mode)
h=(b-a)/N;
t=a:h:b;

if k==2 
  alphap=[0 -1];
  betap=[-0.5 1.5];
  alphac=[0 -1];
  betac=[0 0.5 0.5];
  ep=5/12;ec=-1/12;
  milc=ec/(ep-ec);
elseif k==3 
  alphap=[0 0 -1];
  betap=[5/12 -4/3 23/12];
  alphac=[0 0 -1];
  betac=[0 -1/12 8/12 5/12];
  ep=3/8;ec=-1/24;
  milc=ec/(ep-ec); 
else k==4 
  alphap=[0 0 0 -1];
  betap=[-9/24 37/24 -59/24 55/24];
  alphac=[0 0 0 -1];
  betac=[0 1/24 -5/24 19/24 9/24];
  ep=251/720;ec=-19/720;
  milc=ec/(ep-ec);
end 

if iex==1 
  for j=1:k
     yp(j)=yexact(t(j));
  end 
else
  yp(1:k)=input('Input [y_0 y_1 ... y_(k-1)] ');
end
yp=yp';
yc=yp;
ye=yp;

for i=1:k
  fn(i,1)=f(t(i),yc(i));
  fnex(i,1)=fn(i,1);
  trunc(i)=0;
  mil(i)=0;
end

for n=k:N  
   oldy=yc(n+1-k:n);
   oldf=fn(n+1-k:n); 
   oldyex=ye(n+1-k:n);
   oldfex=fnex(n+1-k:n);
   yp(n+1)=-alphap*oldy+h*(betap*oldf);
   yc(n+1)=yp(n+1);
   for it=1:m
      fc=f(t(n+1),yc(n+1));
      yc(n+1)=-alphac*oldy+h*(betac*[oldf; fc]);
   end  
   mil(n+1)=abs(milc)*abs(yp(n+1)-yc(n+1));
   if mode==1 
     fn(n+1)=f(t(n+1),yc(n+1));
   else
     fn(n+1)=fc;
   end
   ye(n+1)=yexact(t(n+1));
   fnex(n+1)=f(t(n+1),ye(n+1));
   trunc(n+1)=abs(yexact(t(n+1))+alphac*oldyex-h*(betac*[oldfex; fnex(n+1)]));
end

y=yc';