function trapstab(n,el)
% Solves the ivp y'=el y - 1/(1+x)^2 - el/(1+x), y(0)=1 
% by the trapezium rule method
%
fprintf('\n')
%
h=1./n;
h2=h/2.;
hl=abs(h*el);
fprintf('el =%12.6f        h  = %12.6f   hl  = %12.6f\n',el,h,hl)
fprintf('\n')
%
y=1.;
disp('       t          exact         error         relerror')
%
for i=0:n
t=i*h;
t1=(i+1)*h;
ex=(1/(1+t));
er=abs(ex-y);
reler=er/ex;
fprintf('%12.4f%15.6e%15.6e %15.6e\n',t,ex,er,reler)
m= - 1/(1+t)^2 - el/(1+t);
m1=- 1/(1+t1)^2 - el/(1+t1);
a=(1.+el*h2);
b=h2*(m+m1);
c=(1.-el*h2);
y=(y*a+b)/c;
end
