function trap1

%  error as function of M for trap, and both euler's
%  y' = f(t,y)  with   y(0) = y0

global lambda

lambda=10;
y0=0.01;
tmax=1;

clf
% get(gcf)
set(gcf,'Position', [1 1078 573 267])

% exact solution at tmax
a0=(1-y0)/y0;
exact=1/(1+a0*exp(-lambda*tmax));

m=10;
for im=1:6
    m=m*4
    points(im)=m;
    t=linspace(0,tmax,m+1);
    k=t(2)-t(1);
    y_euler=euler(t,y0,k,m+1);
    errorE(im)=abs(exact-y_euler(m+1));
    y_beuler=beuler(t,y0,k,m+1);
    errorBE(im)=abs(exact-y_beuler(m+1));
    y_trap=trap(t,y0,k,m+1);
    errorT(im)=abs(exact-y_trap(m+1));
end

loglog(points,errorE,'--or','MarkerSize',8,'LineWidth',1.3)
hold on
loglog(points,errorBE,'--ok','MarkerSize',8,'LineWidth',1.3)
loglog(points,errorT,'--ob','MarkerSize',8,'LineWidth',1.3)
grid on
xlabel('M')
ylabel('Error')
axis([1 1e5 1e-11 1e-2]);
set(gca,'ytick',[1e-10 1e-8 1e-6 1e-4 1e-2])
set(gca,'YMinorGrid','off')
legend({' Euler',' B Euler',' Trap'},'Location','SouthWest','FontSize',16,'FontWeight','bold')
set(gca,'FontSize',16,'FontWeight','bold')


% euler method
function y=euler(t,y0,k,n)
y=zeros(n,1);
y(1)=y0;
for i=2:n
    y(i)=y(i-1)+k*f(t(i-1),y(i-1));
end

% backward euler method
function ypoints=beuler(t,y0,h,n)
tol=1e-12;
y=y0; fold=f(t(1),y);
ypoints=y0;
for i=2:n
    %  secant method
    yb=y; fb=yb-h*f(t(i),yb)-y;
    yc=y+2*h*(fold+f(t(i),y+h*fold));  fc=yc-h*f(t(i),yc)-y;
    while abs(yc-yb)>tol
        ya=yb; fa=fb;
        yb=yc; fb=fc;
        yc=yb-fb*(yb-ya)/(fb-fa);
        fc=yc-h*f(t(i),yc)-y;
    end
    y=yc; fold=f(t(i),y);
    ypoints=[ypoints, y];
end

% trap method
function ypoints=trap(t,y0,h,n)
tol=1e-12;
y=y0; fold=f(t(1),y);
ypoints=y0;
for i=2:n
    %  secant method
    c=y+0.5*h*fold;
    yb=y;   fb=yb-0.5*h*f(t(i),yb)-c;
    yc=y+0.1*h*fold;   fc=yc-0.5*h*f(t(i),yc)-c;
    while abs(yc-yb)>tol
        ya=yb; fa=fb;
        yb=yc; fb=fc;
        yc=yb-fb*(yb-ya)/(fb-fa);
        fc=yc-0.5*h*f(t(i),yc)-c;
    end
    y=yc; fold=f(t(i),y);
    ypoints=[ypoints, y];
end

% right-hand side of DE
function z=f(t,y)
global lambda
z=lambda*y*(1-y);