alpha = 1.5;  beta = - 2e-1;
t0 = 0;    tf = 30;
y0 = alpha + sin( t0 );
fun = @( t, y ) y * ( cos( t ) / ( alpha + sin( t ) ) + beta );
Nh = 30;
[ tv, uv_forward_euler ] = forward_euler( fun, y0, t0, tf, Nh );
[ tv, uv_heun ] = heun( fun, y0, t0, tf, Nh );
[ tv, uv_runge_kutta_4 ] = runge_kutta_4( fun, y0, t0, tf, Nh );
tv_plot = linspace( t0, tf, 5001 );
y_ex = @( t ) ( alpha + sin( t ) ) .* exp( beta * ( t - t0 ) );
plot( tv_plot, y_ex( tv_plot ), '-k', tv, uv_forward_euler, '.-b', ...
      tv, uv_heun, '.-r', tv, uv_runge_kutta_4, '.-g' );  
grid; axis([-0.1+t0 tf+0.1 -0.25 2.5]);
legend( ' y_{ex}(t)', ' u_{n}, F.E.', ' u_{n}, Heun', ' u_{n}, R.-K. 4' );

errv_n_forward_euler = [];
errv_n_heun = [];
errv_n_runge_kutta_4 = [];
Nhv = [ 30 60 120 240 480 960 ];
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_runge_kutta_4 ] = runge_kutta_4( fun, y0, t0, tf, Nh );    
   % notice that vector index starts from 1 in Matlab
   % while n from zero, n = 0,1,... 
   errv_n_forward_euler = [ errv_n_forward_euler, ...
                       abs(  y_ex( t_bar ) - uv_forward_euler( n + 1 ) ) ];
   errv_n_heun = [ errv_n_heun, ...
                       abs(  y_ex( t_bar ) - uv_heun( n + 1 ) ) ];
   errv_n_runge_kutta_4 = [ errv_n_runge_kutta_4, ...
                       abs(  y_ex( t_bar ) - uv_runge_kutta_4( n + 1 ) ) ];                                                                               
end
% scale factors for visualization of the curves  (h,h), (h,h^2), (h,h^4)
scale_factor_1 = 2 * errv_n_forward_euler( end ) / hv( end );
scale_factor_2 = 2 * errv_n_heun( end ) / hv( end )^2;
scale_factor_3 = 2 * errv_n_runge_kutta_4( end ) / hv( end )^4;
figure;
loglog( hv, errv_n_forward_euler, '-xb', hv, errv_n_heun, '-or', ...
        hv, errv_n_runge_kutta_4, '-dg', ...
        hv, scale_factor_1 * hv, '--k', hv, scale_factor_2 * hv.^2, '.-k',...
        hv, scale_factor_3 * hv.^4, ':k' );
grid; axis( [ 2e-2 2 1e-10 1e0 ])
legend(' e_{n}, F.E.',' u_{n}, Heun', ' e_{n}, R.K.4', ...
       ' h', ' h^2', ' h^4' );
   
conv_ord_Forward_Euler = log( errv_n_forward_euler(end-1) / ...
                              errv_n_forward_euler(end) ) / ...
                                  log( hv(end-1) / hv(end) )
%  conv_ord_Forward_Euler =
%      0.9996

t0 = 0;  tf = 40;
lambda = -0.525;  y0 = 1;
fun = @( t, y ) lambda * y;
Nh = 10; % or 11
[ tv, uv_forward_euler ] = forward_euler( fun, y0, t0, tf, Nh );
[ tv, uv_heun ] = heun( fun, y0, t0, tf, Nh );
[ tv, uv_runge_kutta_4 ] = runge_kutta_4( fun, y0, t0, tf, Nh );
tv_plot = linspace( t0, tf, 5001 );
y_ex = @( t ) y0 * exp( lambda * ( t - t0 ) );
figure;
plot( tv_plot, y_ex( tv_plot ), '-k', tv, uv_forward_euler, '.-b', ...
      tv, uv_heun, '.-r', tv, uv_runge_kutta_4, '.-g' );  
grid; axis([-0.1+t0 tf+0.1 -1.5 2.5 ]);
legend( ' y_{ex}(t)', ' u_{n}, F.E.', ' u_{n}, Heun', ...
     ' u_{n}, R.-K.4',  'Location', 'NorthEast' );
 
h_max_FE_H = 2 / abs( lambda )
%  h_max_FE_H =
%      3.8095
N_h_max_FE_H = ceil( ( tf - t0 ) / h_max_FE_H )
%  N_h_max_FE_H =
%      11
h_max_RK4 = 2.7853 / abs( lambda )
%  h_max_RK4 =
%      5.3053
N_h_max_RK4 = ceil( ( tf - t0 ) / h_max_RK4 )
%  N_h_max_RK4 =
%       8