f = @( x ) x .* log( x ) - sin( x ).^2 ;
df = @( x ) 1 + log( x ) - 2 * sin( x ) .* cos( x );
xnode = 1.9;
df_exact_node = df( xnode )
%  df_exact_node =
%      2.2537
h = 1 / 16;
[ dfh_f_node ] = forward_finite_difference( f, xnode, h )
%  dfh_f_node =
%      2.3178
[ dfh_b_node ] = backward_finite_difference( f, xnode, h )
%  dfh_b_node =
%      2.1861
[ dfh_c_node ] = centered_finite_difference( f, xnode, h )
%  dfh_c_node =
%     2.2519

h_vect = 2.^[ -2 : -1 : -7 ];
err_f = [ ]; err_b = [ ];  err_c = [ ];
for h = h_vect
    [ dfh_f_node ] = forward_finite_difference( f, xnode, h );
    err_f = [ err_f, abs( df_exact_node - dfh_f_node ) ];
    [ dfh_b_node ] = backward_finite_difference( f, xnode, h );
    err_b = [ err_b, abs( df_exact_node - dfh_b_node ) ];
    [ dfh_c_node ] = centered_finite_difference( f, xnode, h );
    err_c = [ err_c, abs( df_exact_node - dfh_c_node ) ];
end
viz_factor = 1.5; % factor to improve the visualization of the line h
loglog( h_vect, err_f, '-ob', h_vect, err_b, '-xr', h_vect, err_c, '-sm', ...
     h_vect, viz_factor * h_vect, '--k', h_vect, h_vect.^2, '-.k' );
grid on;  
xlabel('h'); ylabel('err');
legend( 'e_{+}', 'e_{-}', 'e_{c}', 'h', 'h^2');    
title('Error Finite Differences vs. h, approx df/dx')

conv_order_f = log( err_f( end - 1 ) / err_f( end ) ) / ...
                           log( h_vect( end - 1 ) / h_vect( end ) )
%  conv_order_f =
%       0.9951
conv_order_b = log( err_b( end - 1 ) / err_b( end ) ) / ...
                           log( h_vect( end - 1 ) / h_vect( end ) )
%  conv_order_b =
%       1.0048
conv_order_c = log( err_c( end - 1 ) / err_c( end ) ) / ...
                           log( h_vect( end - 1 ) / h_vect( end ) )
%   conv_order_c =
%       2.0000

a = 3 / 2;  b = 5 / 2;
n = 8; % n is the number of intervals, n + 1 nodes
h = ( b - a ) / n;
x_nodes = linspace( a, b, n + 1 );
x_nodes_internal = x_nodes( 2 : end - 1 );
[ dfh_c_nodes_internal ] = centered_finite_difference( f, x_nodes_internal, h );
dfh_c_node_a = ( - 3 * f( a ) + 4 * f( a + h ) - f( a + 2 * h ) ) / ( 2 * h );
dfh_c_node_b = ( 3 * f( b ) - 4 * f( b - h ) + f( b - 2 * h ) ) / ( 2 * h );
dfh_c_nodes = [ dfh_c_node_a, dfh_c_nodes_internal, dfh_c_node_b ];
x_values = linspace( a, b , 100 * ( n + 1 ) ); % for visualization of exact values
df_exact_values = df( x_values );
figure; plot( x_nodes, dfh_c_nodes, 'xr', x_values, df_exact_values, '-k' );
legend( '(d_c f)(x_i)', 'df/dx(x)');