function AB4(n)
% Solves the ivp y'= -y - 1/(1+x)^2 + 1/(1+x), y(0)=1 
% by the four-step method A-B method
% Y_{n+1} - Y_n  =  (h/24) (55 Y_{n}' -59  Y_{n-1}' + 37 Y_{n-2}' - 9 Y_{n-3}')
%
fprintf('\n')
%
for iit=1:10
%
nn=2^(iit-1)*n;
uu(iit)=1/nn;
%
h=1./nn;
%
h24=h/24;
%
y0=1.;
y1=1/(1+h);
y2=1/(1+2*h);
y3=1/(1+3*h);
%
er=zeros(nn,1);
for i=5:nn
t=(i-1)*h;
t1=(i-2)*h;
t2=(i-3)*h;
t3=(i-4)*h;
t4=(i-5)*h;
ex=1/(1+t);
m1= -y3- 1/(1+t1)^2 + 1/(1+t1);
m2= -y2- 1/(1+t2)^2 + 1/(1+t2);
m3= -y1- 1/(1+t3)^2 + 1/(1+t3);
m4= -y0 - 1/(1+t4)^2 + 1/(1+t4);
y4=y3+h24*(55*m1-59*m2+37*m3-9*m4);
er(i)=abs(ex-y4);
y0=y1;
y1=y2;
y2=y3;
y3=y4;
end
u1(iit)=max(er);
fprintf('h =%12.6f        maxerr = %15.6e\n',h,u1(iit))
end
fprintf('\n')
for iit=1:9
      p=log(u1(iit)/u1(iit+1))/log(2.);
      fprintf('h =%12.6f        rate %15.6e\n',uu(iit),p)
end
loglog(uu,u1,'b+')
xlabel('h');ylabel('Error');
