function logistic_error

%  max error for various DE solvers of logistic equation
%  y' = ry(1-y)  with   y(0) = y0

%  suggestion:  in trap function (see below) first use  tol=1.0e-7  then use  tol=1.0e-8

% clear all previous variables and plots
clear *
clf


% set parameters for problem
t0=0; y0=0.1;
tmax=1;


% loop that increases M 
for im=1:4

	M=10^im;
	points(im)=M;
	t=linspace(t0,tmax,M+1);
	h=t(2)-t(1);


	y_exact=exact(t,y0,M+1);

	% euler
        y_euler=euler('de_f',t,y0,h,M+1);
	e_euler(im)=norm(y_exact-y_euler,inf);


	% backward euler
        y_beuler=beuler('de_f',t,y0,h,M+1);
	e_beuler(im)=norm(y_exact-y_beuler,inf);
   	
	%  trap
        y_trap=trap('de_f',t,y0,h,M+1);
        e_trap(im)=norm(y_exact-y_trap,inf);

    	%  RK4
        y_rk4=rk4('de_f',t,y0,h,M+1);
        e_rk4(im)=norm(y_exact-y_rk4,inf);

end;

% plot error curves
loglog(points,e_euler,'--ro','MarkerSize',10)
hold on
loglog(points,e_beuler,'--rs','MarkerSize',10)
grid on
loglog(points,e_trap,'--m*','MarkerSize',10)
loglog(points,e_rk4,'--bd','MarkerSize',10)

% define legend and axes used in plot
legend(' Euler',' B Euler',' Trap',' RK4',3);
xlabel('Number of Grid Points Used to Reach t=1','FontSize',14,'FontWeight','bold')
ylabel('Error','FontSize',14,'FontWeight','bold')
title('Logistic Equation: Max Error','FontSize',14,'FontWeight','bold')
set(gca,'ytick',[1e-15  1e-13 1e-11 1e-9 1e-7 1e-5 1e-3 1e-1]);

% have MATLAB use certain plot options (all are optional)
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  

hold off


% right-hand side of DE
function z=de_f(t,y)
r=10;
z=r*y*(1-y);

% exact solution

function y=exact(t,y0,n)
a0=(1-y0)/y0;
r=10;
y=y0;
for i=2:n
	yy=1/(1+a0*exp(-r*t(i)));
	y=[y, yy];
end;


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


% trap method
function ypoints=trap(f,t,y0,h,n)
%  suggestion:  first use  tol=1.0e-7  then use  tol=1.0e-8

tol=1.0e-7;
y=y0; fold=feval(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*feval(f,t(i),yb)-c;
	yc=y+0.1*h*fold;  fc=yc-0.5*h*feval(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*feval(f,t(i),yc)-c;
	end;
	y=yc; fold=feval(f,t(i),y);
	ypoints=[ypoints, y];
end;


% euler method
function ypoints=euler(f,t,y0,h,n)
y=y0;
ypoints=y0;
for i=2:n
	yy=y+h*feval(f,t(i-1),y);
	ypoints=[ypoints, yy];
	y=yy;
end;

% backward euler method

function ypoints=beuler(f,t,y0,h,n)
tol=1.0e-7;
y=y0; fold=feval(f,t(1),y);
ypoints=y0;
for i=2:n
	%  secant method
	yb=y;         fb=yb-h*feval(f,t(i),yb)-y;
	yc=y+2*h*(fold+feval(f,t(i),y+h*fold));  fc=yc-h*feval(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*feval(f,t(i),yc)-y;
	end;
	y=yc; fold=feval(f,t(i),y);
	ypoints=[ypoints, y];
end;