alpha = pi/2;  beta = pi/3;
fun = @( t, y ) alpha * y * ( 1 - 1 / beta * y );
dfun_y = @( t, y ) alpha * ( 1 - 2 / beta * y );
t0 = 0;  tf = 20;  y0 = 0.4;
Nh = 20;  
[ tv, uv_forward_euler ] = forward_euler( fun, y0, t0, tf, Nh );
[ tv, uv_backward_euler ] = backward_euler( fun, dfun_y, y0, t0, tf, Nh );
tv_plot = linspace( t0, tf, 1001 );
y_ex = @( t ) beta * exp( alpha * ( t - t0 ) + log( y0 / ( beta - y0 ) ) ) ./  ...
           ( 1 + exp( alpha * ( t - t0 ) + log( y0 / ( beta - y0 ) ) ) );
figure
plot( tv_plot, y_ex( tv_plot ), '-k', ...
      tv, uv_forward_euler, '-xb', ...
      tv, uv_backward_euler, '-sr' );
grid on;  xlabel('t');  ylabel('y');
legend(' y_{ex}(t)',' u_{n}, F.E.', ' u_{n}, B.E.' );

y = linspace( y0, beta, 1001 );
df_val = alpha * ( 1 - 2 / beta * y );
df_max = max( abs( df_val ) )
% df_max =
%     1.5708
h_max = 2 / df_max
% h_max =
%     1.2732
Nh_max = ceil( ( tf - t0 ) / h_max )
% Nh_max =
%     16
       
Nh = 15;  % unstable 
[ tv, uv_forward_euler ] = forward_euler( fun, y0, t0, tf, Nh );
[ tv, uv_backward_euler ] = backward_euler( fun, dfun_y, y0, t0, tf, Nh );
tv_plot = linspace( t0, tf, 1001 );
y_ex = @( t ) beta * exp( alpha * ( t - t0 ) + log( y0 / ( beta - y0 ) ) ) ./  ...
              ( 1 + exp( alpha * ( t - t0 ) + log( y0 / ( beta - y0 ) ) ) );
figure
plot( tv_plot, y_ex( tv_plot ), '-k', ...
      tv, uv_forward_euler, '-xb', ...
      tv, uv_backward_euler, '-sr' );  
grid on;  xlabel('t');  ylabel('y');
legend(' y_{ex}(t)',' u_{n}, F.E.', ' u_{n}, B.E.' );