f = @(x) exp(-x.^2 / 2);
a = -5; b = 5;
x_values = linspace( a, b, 1001 );
f_values = f( x_values );
n_vect = [2 7 12 22 27 32 ];
for n = n_vect
   x_nodes = linspace( a, b, n + 1 );
   y_nodes = f( x_nodes );
   yp1_values = interp1( x_nodes, y_nodes, x_values ); 
   plot( x_values, f_values, '--k', x_values, yp1_values, '-k', ...
           x_nodes, y_nodes, 'xk');
   legend( 'f(x)', '\Pi_1^H f(x)', '(x_i,f(x_i))');    
   pause(0.1)    
end

n_vect = 2:32;
err_yp1 = [];
for n = n_vect
   x_nodes = linspace( a, b, n + 1 );
   y_nodes = f( x_nodes );
   yp1_values = interp1( x_nodes, y_nodes, x_values ); 
   err_yp1 = [ err_yp1, max( abs( f_values - yp1_values ) ) ];
end
figure; semilogy( n_vect, err_yp1, '-ok' );
grid on

n = 7;
x_nodes = linspace( a, b, n + 1 );
y_nodes = f( x_nodes );
yp1_values = interp1( x_nodes, y_nodes, x_values ); 
s3_values = spline( x_nodes, y_nodes, x_values );
figure; 
plot( x_values, f_values, '--k', x_values, yp1_values, '-k', ...
        x_values, s3_values, '-r', x_nodes, y_nodes, 'xk' );
legend( 'f(x)', '\Pi_1^H f(x)', 's_3(x)', '(x_i,f(x_i))' );    

n_vect = 2:32;
err_s3 = [];
for n = n_vect
   x_nodes = linspace( a, b, n + 1 );
   y_nodes = f( x_nodes );
   s3_values = spline( x_nodes, y_nodes, x_values ); 
   err_s3 = [ err_s3, max( abs( f_values - s3_values ) ) ];
end
figure; semilogy( n_vect, err_s3, '-ok' );
grid on