function beuler0

%  backward euler method for various M values
%  compares with exact solution
%  uses either exact solution of difference equation or secant solutiuon
%  y' = f(t,y)  with   y(0) = y0

global lambda

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

%  icase = 1: use exact solution of difference equation
%  icase = 2: use secant method solution of difference equation
icase=1;

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

% calculate and plot exact solution
tt=linspace(0,tmax,100);
a0=(1-y0)/y0;
for it=1:100
    exact(it)=1/(1+a0*exp(-lambda*tt(it)));
end
plot(tt,exact,'r','LineWidth',1.8)
hold on
grid on
box on
xlabel('t-axis')
ylabel('Solution')
set(gca,'FontSize',16,'FontWeight','bold')
axis([0 1 0 1.05])

% compute numerical solution and plot
m=1;
for im=1:3
    m=4*m
    t=linspace(0,tmax,m+1);
    k=t(2)-t(1);
    
    % backward Euler method
    if icase==1
        y_beuler=bbeuler(t,y0,k,m+1);
    else
        y_beuler=beuler(t,y0,k,m+1);
    end
    
    if im==1
        plot(t,y_beuler,'--ks','MarkerSize',9,'LineWidth',1.3)
        legend({' Exact',' M = 4'},'Location','NorthWest','FontSize',16,'FontWeight','bold')
        pause
    elseif im==2
        plot(t,y_beuler,'--m*','MarkerSize',9,'LineWidth',1.3)
        legend({' Exact',' M = 4',' M = 16'},'Location','NorthWest','FontSize',16,'FontWeight','bold')
        pause
    else
        plot(t,y_beuler,'--bo','MarkerSize',9,'LineWidth',1.3)
        legend({' Exact',' M = 4',' M = 16',' M = 64 '},'Location','NorthWest','FontSize',16,'FontWeight','bold')
    end
    
    
end

% backward euler method using exact soltion of difference equation
function y=bbeuler(t,y0,k,n)
global lambda
y=zeros(n,1);
kl=k*lambda;
y(1)=y0;
for i=2:n
    y(i)=(-1+kl+sqrt((1-kl)^2+4*kl*y(i-1)))/(2*kl);
end

% backward euler method
function ypoints=beuler(t,y0,h,n)
tol=1e-6;
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

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