function multi1(n)
% Solves the ivp y'= -y - 1/(1+x)^2 + 1/(1+x), y(0)=1 
% by the explicit three-step method
% Y_{n+1} + 9 Y_n - 9 Y_{n-1} - Y_{n-2} = 6 h (Y_n' + Y_{n-1}')
%
fprintf('\n')
%
for iit=1:6
%
nn=2^(iit-1)*n;
uu(iit)=1/nn;
%
h=1./nn;
%
y0=1.;
y1=1./(1.0+h);
y2=1./(1.0+2.0*h);
%
for i=2:nn
t=i*h;
t1=(i-1)*h;
ex=(1/(1+t));
er(i)=abs(ex-y2);
m= -y2- 1/(1+t)^2 + 1/(1+t);
m1= -y1 - 1/(1+t1)^2 + 1/(1+t1);
y3=-9.0*y2+9.0*y1+y0+6.0*h*(m+m1);
y0=y1;
y1=y2;
y2=y3;
%fprintf('point =%12.6f        error = %15.6e\n',t,er(i))
end
u1(iit)=max(er);
fprintf('h =%12.6f        maxerr = %15.6e\n',h,u1(iit))
end
fprintf('\n')
for iit=1:5
      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');
