t0 = 0;   tf = 5;
y0 = [ 1; 0.6 ];
A = - [ 3 1; 1 1 ];
g = @( t ) [ 0; 0 ] * t;
fun = @( t, y ) A * y;
Nh = 25; % 8 % 9

[ tv, uv_forward_euler ] = forward_euler_system( fun, y0, t0, tf, Nh );
figure;
plot( tv, uv_forward_euler( 1, : ), '.-b', ...
          tv, uv_forward_euler( 2, : ), '.-r' );
grid; axis([-0.1+t0 tf+0.1 -1 1.25]);
title('System of ODEs, Forward Euler');
legend('y(1)','y(2)','Location','NorthEast');

[ tv, uv_backward_euler ] = backward_euler_system_nhcc( A, g, y0, t0, tf, Nh );
figure;
plot( tv, uv_backward_euler( 1, : ), '.-b', ...
          tv, uv_backward_euler( 2, : ), '.-r' );
grid;  axis([-0.1+t0 tf+0.1 -1 1.25]);
title('System of ODEs, Backward Euler');
legend('y(1)','y(2)','Location','NorthEast');

lambda = ( eig( A ) )'
%  lambda =
%     -3.4142    -0.5858
h_max = 2 / max( abs( lambda ) )
%  h_max =
%     0.5858
N_h_max = ceil( ( tf - t0 ) / h_max )
%  N_h_max =
%     9