function ivp2
 
%  comparison of various DE solvers
%  y' = f(t,y)  with   y(0) = y0             '
 
% clear all previous variables and plots 
clear *
clf
 
% set parameters for problem
M=3200;
y0=2;
tmax=3;
 
%  calculate time points
t=linspace(0,tmax,M+1);
h=t(2)-t(1);
 
fprintf('\n    backward euler, trapezoidal, and RK4 with M = %3.0f \n\n',M)    
 
% backward euler
y_beuler=beuler('de_f',t,y0,h,M+1);
plot(t,y_beuler,'--bo')
hold on
 
%  trap
y_trap=trap('de_f',t,y0,h,M+1);
plot(t,y_trap,'--rs')
 
%  RK4
y_rk4=rk4('de_f',t,y0,h,M+1);
plot(t,y_rk4,'--b*')
 
legend(' B Euler', ' Trap', ' RK4',1);
title(['M = ',num2str(M)],'FontSize',14,'FontWeight','bold')
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
 
% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 12 for the plot
set(gca,'FontSize',12);
% 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)
% exact sol  y=cos(w*t)*(1+exp(-b*t))
b=40; w=3;
z=y*(-w*tan(w*t)-b/(1+exp(b*t)));
 
 
% 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)
tol=0.0000001;
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;
 
% backward euler method
function ypoints=beuler(f,t,y0,h,n)
tol=0.0000001;
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;