function ivp3

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0        '
%  where y and f are scalar functions

% clear all previous variables and plots
clear *
clf

% time interval
tmax=3;

%  initial value
y0=2;

% calculate, and plot, exact solution
tt=linspace(0,tmax,200);
b=40; w=3;
for it=1:200
	exact(it)=cos(w*tt(it))*(1+exp(-b*tt(it)));
end;
plot(tt,exact,'k')
pause
hold on

%  calculate solution using a MATLAB routine
%[t,y] = ode45(@rhs,[0 tmax],y0);  titleX=' ode45';
[t,y] = ode23s(@rhs,[0 tmax],y0);  titleX=' ode23s';

%  plot solution as function of t
plot(t,y,'r')
legend(' Exact',titleX,4);

% commands to label each axes
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('y-axis','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
box on
% 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'); 

% commands used to automatically size window
% get(gcf)
% set(gcf,'Position', [475 511 484 278]);

hold off


%  define f(t,y)
function dy=rhs(t,y)
% exact sol  y=cos(w*t)*(1+exp(-b*t))
b=40; w=3;
f=-w*tan(w*t)-b/(1+exp(b*t));
dy=f*y;