function [ ] = prob_2_19
% 
% 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 epsilon kappa;  
epsilon = 1/98;
kappa = 3;

options = odeset('Stats','on');

y0 = [ 1, 0];

disp( 'Solving with ode15s' );
[ts,ys] = ode15s( @f, [0, 240], y0, options );
figure;
semilogx( ts, ys(:,2), '-og' );
title( 'ODE15S solution' );
axis( [ 0.01 100 0 1 ] );

disp( 'Solving with ode45' );
[ts,ys] = ode45( @f, [0, 240], y0, options );
figure;
semilogx( ts, ys(:,2), '-og' );
title( 'ODE45 solution' );
axis( [ 0.01 100 0 1 ] );

return;



function [dydt] = f(t,y)
%
% y is 2x1
%
global epsilon kappa;  
  
dydt =             [ -y(1) - y(1)*y(2) + epsilon*kappa*y(2);
       (1/epsilon)*(  y(1) - y(1)*y(2) - epsilon*kappa*y(2) ) ];

return;