function rk4

%  RK4
%  y' = f(t,y)  with   y(0) = y0

global lambda

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


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

% RK4 method
function ypoints=rk4s(t,y0,h,n)
y=y0;
ypoints=y0;
for i=2:n
    k1=h*f(t(i-1),y);
    k2=h*f(t(i-1)+0.5*h,y+0.5*k1);
    k3=h*f(t(i-1)+0.5*h,y+0.5*k2);
    k4=h*f(t(i),y+k3);
    yy=y+(k1+2*k2+2*k3+k4)/6;
    ypoints=[ypoints, yy];
    y=yy;
end