function eulers

%  euler method for various M values when solving
%  y' = ry(1-y)  with   y(0) = y0      '

% clear all previous variables and plots
clear *
clf

% parameters for calculation
global r;
r=50 ;
t0=0; y0=0.01;
tmax=1;

% set the position and size of the plot
% get(gcf)
% set(gcf,'Position', [1203 732 515 307]);

% calculate and plot exact solution
tt=linspace(t0,tmax,100);
a0=(1-y0)/y0;
for it=1:100
	exact(it)=1/(1+a0*exp(-r*tt(it)));
end;

% loop used to increase M value
M=1;
for icounter=1:4
	M=M*4

	% calculate Euler solution
	t=linspace(t0,tmax,M+1);
	h=t(2)-t(1);
	euler=y0;
	y=y0;
	for i=2:M+1
		yy=y+r*h*y*(1-y);
		euler=[euler, yy];
		y=yy;
	end;

	% calculate backward euler solution
	y_beuler=beuler('de_f',t,y0,h,M+1);

	% plot solutions
	clf
	hold on
	plot(tt,exact,'k','LineWidth',1)
	plot(t,euler,'--ro','MarkerSize',7)
	plot(t,y_beuler,'--bo','MarkerSize',7)
	legend(' Exact',' Euler',' B Euler',2)
	set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold')
	title(['M = ',num2str(M)],'FontSize',14,'FontWeight','bold')

	% define axes and title used in plot
	xlabel('t-axis','FontSize',14,'FontWeight','bold')
	ylabel('Solution','FontSize',14,'FontWeight','bold')

	% have MATLAB use certain plot options (all are optional)
	box on
%	axis([0 1 0 1.05])
%	set(gca,'ytick',[0 0.2 0.4 0.6 0.8 1.0]);
	set(gca,'FontSize',14); 

	pause
	hold off
end
 


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

% backward euler method
function ypoints=beuler(f,t,y0,h,n)
global r;
tol=0.00001;
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;
	yc=(-1+r*h+sqrt((1-r*h)^2+4*r*h*y))/(2*r*h);
	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;