C1 = 0.15;    C2 = 0.075;
b1 = 0.002;   b2 = 0;
d1 = 0.0210;  d2 = 0.0325;
t0 = 0;       tf = 600;
y0 = [ 55; 20 ]; % [35; 40];
fun = @( t, y ) [ C1 * y( 1 ) * ( 1 - b1 * y( 1 ) - d2 * y( 2 ) ); ...
                     - C2 * y( 2 ) * ( 1 - b2 * y( 2 ) - d1 * y( 1 ) ) ];
Nh = 5000; 
[ tv, uv_forward_euler ] = forward_euler_system( fun, y0, t0, tf, Nh );
plot( tv, uv_forward_euler( 1, : ), '.-b', tv, uv_forward_euler( 2, : ), '.-r' );
grid; axis([-0.1+t0 tf+0.1 -1 80]);
legend('Prey','Predator','Location','SouthEast');

equilibrium_point = [ 47.6190; 27.8388 ];
figure; plot( uv_forward_euler( 1, : ), uv_forward_euler( 2, : ), '-k', ...
      equilibrium_point( 1 ), equilibrium_point( 2 ), '+r' );
xlabel('Prey'); ylabel('Predator');
grid; axis([10 80 10 50]);