function euler0

%  euler method for various M values
%  compares with exact solution
%  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])

% 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);
    
    % use Euler method
    y_euler=euler(t,y0,k,m+1);
    
    if im==1
        plot(t,y_euler,'--ks','MarkerSize',9,'LineWidth',1.3)
        legend({' Exact',' M = 4'},'Location','NorthWest','FontSize',16,'FontWeight','bold')
        pause
    elseif im==2
        plot(t,y_euler,'--m*','MarkerSize',9,'LineWidth',1.3)
        legend({' Exact',' M = 4',' M = 16'},'Location','NorthWest','FontSize',16,'FontWeight','bold')
        pause
    else
        plot(t,y_euler,'--bo','MarkerSize',9,'LineWidth',1.3)
        legend({' Exact',' M = 4',' M = 16',' M = 64 '},'Location','NorthWest','FontSize',16,'FontWeight','bold')
    end
    
    
end


% 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


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