%
%  hw04.m -- problem 9.1, heath, page 297.
%
%  2nd-order runge-kutta
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  set parameters, IVs & coefficients

%N = 1000;  L = 10;  dt = L/N;
%y0 = 0.5;  z0 = 1.0;

N = 200;  L = 10;  dt = L/N;
y0 = 5.0;  z0 = 1.0;

b  = 1.0;  c  = 1.0;  d = 10.0;

t  = 0:dt:L;
ys = zeros(1,N+1);
zs = zeros(1,N+1);
check = zeros(1,N+1);

%  initial condition & fwd euler step (change 0 to 1 for no osc'ns)

ys(1) = y0;  zs(1) = z0;
check(1) = c*ys(1) - d*log(ys(1)) + c*zs(1) - b*log(zs(1));
disp(['first check = ' num2str(check(1))])

%  iteration loop

for j=2:1:N+1
  ky1 = dt*( b*(ys(j-1)) - c*ys(j-1)*zs(j-1));
  kz1 = dt*(-d*(zs(j-1)) + c*ys(j-1)*zs(j-1));

  ky2 = dt*( b*(ys(j-1) + ky1) - c*(ys(j-1) + ky1)*(zs(j-1) + kz1));
  kz2 = dt*(-d*(zs(j-1) + kz1) + c*(ys(j-1) + ky1)*(zs(j-1) + kz1));

  ys(j) = ys(j-1) + (ky1 + ky2)/2;
  zs(j) = zs(j-1) + (kz1 + kz2)/2;

  check(j) = c*ys(j) - d*log(ys(j)) + c*zs(j) - b*log(zs(j));
end 

disp(['last check  = ' num2str(check(N+1))])

max(check)-min(check)

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
%  solution plot

clf;subplot(2,1,1)
plot(t,ys,'k.')
hold on
plot(t,ys,'b')
plot(t,zs,'k.')
plot(t,zs,'g')
title('\bf predator/prey model ODEs')
xlabel('\bf t-axis');  ylabel('\bf y,z-axis');

subplot(2,1,2)
plot(t,check,'k.')
title('\bf predator/prey model ODEs')
xlabel('\bf t-axis');  ylabel('\bf check-axis');