fun = @( t, y ) 1 - y.^2;
t0 = 0;   tf = 5;
y0 = ( exp( 1 ) - 1 ) / ( exp( 1 ) + 1 );
Nh = 10;
[ tv, uv_forward_euler ] = forward_euler( fun, y0, t0, tf, Nh );
[ tv, uv_heun ] = heun( fun, y0, t0, tf, Nh );
tv_plot = linspace( t0, tf, 1001 );
y_ex = @( t ) ( exp( 2 * t + 1 ) - 1 ) ./ ( exp( 2 * t + 1 ) + 1 );
figure
plot( tv_plot, y_ex( tv_plot ), '-k', ...
      tv, uv_forward_euler, '-xb', tv, uv_heun, '-or' );
grid on; xlabel('t'); ylabel('y');
legend(' y_{ex}(t)',' u_{n}, F.E.',' u_{n}, Heun' );

lambda = - 0.5;
fun = @( t, y ) lambda * y;
t0 = 0;  tf = 15;
y0 = 1;
Nh = 10;
[ tv, uv_forward_euler ] = forward_euler( fun, y0, t0, tf, Nh );
[ tv, uv_heun ] = heun( fun, y0, t0, tf, Nh );
tv_plot = linspace( t0, tf, 1001 );
y_ex = @( t ) y0 * exp( lambda * t );
figure
plot( tv_plot, y_ex( tv_plot ), '-k', ...
      tv, uv_forward_euler, '-xb', tv, uv_heun, '-or' );
grid on
xlabel('t');
ylabel('y');
legend(' y_{ex}(t)',' u_{n}, F.E.',' u_{n}, Heun' );

[ tv, uv_backward_euler ] = backward_euler_modelproblem( lambda, ...
                                                          y0, t0, tf, Nh );
[ tv, uv_crank_nicolson ] = crank_nicolson_modelproblem( lambda, ...
                                                          y0, t0, tf, Nh );
figure
plot( tv_plot, y_ex( tv_plot ), '-k', ...
      tv, uv_backward_euler, '-dg', tv, uv_crank_nicolson, '-sm'  );
grid on; xlabel('t'); ylabel('y');
legend(' y_{ex}(t)',' u_{n}, B.E.',' u_{n}, C.-N.' );

errv_n_forward_euler = [];
errv_n_backward_euler = [];
errv_n_heun = [];
errv_n_crank_nicolson = [];
Nhv = [ 15 30 60 120 240 480 ];
hv = ( tf - t0 ) ./ Nhv;
t_bar = 10;
for Nh = Nhv
    h = ( tf - t0 ) / Nh;   % h
    n = ( t_bar - t0 ) / h; % step n varies with Nh
   [ tv, uv_forward_euler ] = forward_euler( fun, y0, t0, tf, Nh );
   [ tv, uv_heun ] = heun( fun, y0, t0, tf, Nh );
   [ tv, uv_backward_euler ] = backward_euler_modelproblem( lambda, ...
                                                         y0, t0, tf, Nh );
   [ tv, uv_crank_nicolson ] = crank_nicolson_modelproblem( lambda, ...
                                                         y0, t0, tf, Nh );                                                      
   errv_n_forward_euler = [ errv_n_forward_euler, ...
                       abs( uv_forward_euler( n + 1 ) - y_ex( t_bar ) ) ];
   errv_n_heun = [ errv_n_heun, ...
                       abs( uv_heun( n + 1 ) - y_ex( t_bar ) ) ];
   errv_n_backward_euler = [ errv_n_backward_euler, ...
                       abs( uv_backward_euler( n + 1 ) - y_ex( t_bar ) ) ];
   errv_n_crank_nicolson = [ errv_n_crank_nicolson, ...
                       abs( uv_crank_nicolson( n + 1 ) - y_ex( t_bar ) ) ];                                                                               
end
figure
loglog( hv, errv_n_forward_euler, '-xb', ...
        hv, errv_n_heun, '-or', ...
        hv, errv_n_backward_euler, '-dg', ...
        hv, errv_n_crank_nicolson, '-sm', ...
        hv, 1.5 * errv_n_backward_euler( 1 ) / hv( 1 ) * hv, '--k',...
        hv, 1.5 * errv_n_heun( 1 ) / hv( 1 )^2 * hv.^2, '--k' );
grid on; xlabel('h'); ylabel('e_n');
legend(' e_{n}, F.E.',' u_{n}, Heun', ' e_{n}, B.E.',' e_{n}, C.-N.', ...
       ' h', ' h^2' );