function [ ] = prob_1_15
% 
% 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.
% 
%-----

% Parameters:
global a c;
a=2; 
c=1;
% Initial conditions:
t(1) = 0; x(1) = 1; y(1) = 3; 
% Final time:
tmax = 10.0;

% A sample of step sizes:
hArray = 10.^[-1:-1:-4];
%hArray = 10.^-4;

for hi=hArray,
  disp( [ 'hi = ', num2str(hi) ] );
  t(2:end) = []; x(2:end) = []; y(2:end) = [];
  [ t, x, y ] = eulerMethod( hi, t, x, y, tmax );
  figure; plot( x, y, '-o' );
  plotBigG( t, x, y, hi );
end

return;


%------------------------------------------------------------------------
%------------------------------------------------------------------------

function [ t, x, y ] = eulerMethod( hi, t, x, y, tmax )

global a c;

%while( max(t) < tmax )
while( t(end) < tmax )
  t(end+1) = t(end) + hi;
  %x(end+1) = x(end) + hi*a*x(end)*(1-y(end));
  %y(end+1) = y(end) - hi*c*x(end)*(1-y(end));
  x(end+1) = x(end) + hi*a*x(end)*(1-y(end));
  y(end+1) = y(end) - hi*c*y(end)*(1-x(end));
end

return;


%------------------------------------------------------------------------
%------------------------------------------------------------------------

function [] = plotBigG( t, x, y, hi )

global a c;

G = x.^(-c) .* y.^(-a) .* exp( c*x + a*y );
figure; plot( t, G, '-o' );
xlabel( 'time (s)' ); ylabel( 'G (unitless)' );
title( [ 'hi = ', num2str(hi) ] );

return;