function fig7_6

%  comparison of trap, euler and b euler
%  y' = ry(1-y)  with   y(0) = y0 

M=30;
r=10;
t0=0; y0=0.01;
tmax=3;
t=linspace(t0,tmax,M+1);

% 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

% euler
h=t(2)-t(1);
hr=h*r;
y=y0;
euler=y0;
for i=2:M+1
    yy=y+hr*y*(1-y);
    euler=[euler, yy];
    y=yy;
end

% backward euler
y=y0;
beuler=y0;
for i=2:M+1
    yy=-0.5*(1-hr-sqrt((1-hr)^2+4*hr*y))/hr;
    beuler=[beuler, yy];
    y=yy;
end

% trapezoidal
y=y0;
trap=y0;
for i=2:M+1
    b=2*hr*y*(1+0.5*hr*(1-y));
    yy=(-1+0.5*hr+sqrt((1-0.5*hr)^2+b))/hr;
    trap=[trap, yy];
    y=yy;
end


% plot results
clf
% get(gcf)
set(gcf,'Position', [9 870 742 475])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

subaxis(2,1,1,1,'MT',-0.01,'MB',0.075,'MR',-0.02,'ML',0.05,'P',0.03)
hold on
box on
grid on
axis([0 tmax 0 1.1])
plot(tt,exact,'k','LineWidth',1.2)
plot(t,euler,'--bo','MarkerSize',8,'LineWidth',1.6)
plot(t,beuler,'--d','color', [0 0.5 0],'MarkerSize',8,'LineWidth',1.6)
plot(t,trap,'--r*','MarkerSize',8,'LineWidth',1.6)
set(gca,'ytick',[0 0.5 1])
set(gca,'xtick',[0  1  2 3])

legend({' Exact',' Euler',' B Euler',' Trap'},'Location','SouthEast','FontSize',16,'FontWeight','bold')

%xlabel('t-axis')
ylabel('Solution')
set(gca,'FontSize',16,'FontWeight','bold')
hold off

subaxis(2,1,1,2)
M=30;
r=10;
t=linspace(0,tmax,M+1);
% leapfrog solution
h=t(2)-t(1);
hr=h*r;
y=y0;
yy=1/(1+a0*exp(-hr));
leap=[y, yy];
for i=3:M+1
	yyy=y+2*hr*yy*(1-yy);
	leap=[leap, yyy];
	y=yy;
	yy=yyy;
end

plot(tt,exact,'k','LineWidth',1.2)
hold on
grid on
plot(t,leap,'--*b','LineWidth',1.6)

M=150;
r=10;
t=linspace(0,tmax,M+1);
% leapfrog solution
h=t(2)-t(1);
hr=h*r;
y=y0;
yy=1/(1+a0*exp(-hr));
leap=[y, yy];
for i=3:M+1
	yyy=y+2*hr*yy*(1-yy);
	leap=[leap, yyy];
	y=yy;
	yy=yyy;
end
plot(t,leap,'-r','LineWidth',1.6)

legend({' Exact',' Leap (M=30)',' Leap (M=150)'},'Location','NorthWest','FontSize',16,'FontWeight','bold')
set(gca,'ytick',[-0.5 0 0.5 1 1.5 2])
axis([0 tmax -0.5 2])
set(gca,'xtick',[0  1  2 3])
xlabel('t-axis')
ylabel('Solution')
set(gca,'FontSize',16,'FontWeight','bold')                       		
hold off

%exportgraphics(gcf,'/Users/holmes/Desktop/etal.eps')


















