function [ ] = prob_2_22
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-07
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

global f_2 f3 omega_2 omega_3 zeta;
f_2 = 0.05;
omega_2 = 1;
omega_3 = 2;
zeta = 0.1;

options = odeset('RelTol',1e-8,'AbsTol',1e-8);

y0 = [ 0; -0.6  ]; f3 = 0.5;
solveAndPlot( y0, f3, options );

y0 = [ 0; -0.8  ]; f3 = 0.92;
solveAndPlot( y0, f3, options );

y0 = [ 0; -0.59 ]; f3 = 0.99;
solveAndPlot( y0, f3, options );

y0 = [ 0; -0.8  ]; f3 = 1.005;
solveAndPlot( y0, f3, options );

y0 = [ 0;  0.0  ]; f3 = 1.35;
solveAndPlot( y0, f3, options );

return;



function [] = solveAndPlot( y0, f3, options )

[ts,ys] = ode45( @f, [0, 500], y0, options );
ndx = find( ts >= 400 );
figure;
plot( ys(ndx,1), ys(ndx,2), '-g' );
title( [ 'y0(1) = ', num2str(y0(1)), ' y0(2) = ', num2str(y0(2)), ...
	 ' f3 = ', num2str(f3) ] );
drawnow;

return;



function [dydt] = f(t,y)
%
% y is 2x1
%

global f_2 f3 omega_2 omega_3 zeta;
  
dydt = [ 
          y(2),
         -2*zeta*y(2) - y(1) - 1.5*y(1)^2 - 0.5*y(1)^3 + ...
        f_2*cos(omega_2*t) + f3*cos(omega_3*t);
       ];

return;