function [ ] = prob_2_30
% 
% 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('Mass',@mass,...
		 'MassSingular','no','MStateDependence','none',...
		 'Stats','on');

y0 = [ 1, 0];

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

% I don't know why this ode solver has so much problems with this ODE?
[ts,ys] = ode15s( @f, [0, 240], y0, options );
figure;
semilogx( ts, ys(:,2), '-og' );
axis( [ 0.01 100 0 1 ] );

return;



function [A] = mass(t,y)
  
global epsilon kappa;  

A = diag(1,epsilon);

return;



function [dydt] = f(t,y)
%
% y is 2x1
%
global epsilon kappa;

dydt = [ -y(1) - y(1)*y(2) + epsilon*kappa*y(2);
       (  y(1) - y(1)*y(2) - epsilon*kappa*y(2) ) ];

return;