A = [ 4 -2 -1; -1 3 -1; -1 -3 5 ];
x_ex = ones( 3, 1 );
b = A * x_ex;
% Note: MATLAB may perfom pivoting even when not strictly necessary
[ L, U, P ] = lu( A );
[ y ] = forward_substitutions( L, P * b );
[ x ] = backward_substitutions( U, y );
err = norm( x - x_ex )
%  err =
%    3.8459e-16

P
% P =
%    1     0     0
%    0     0     1
%    0     1     0

n = 20;
A = diag( 4 * ones( n, 1 ), 0 ) + ...
    diag( -1 * ones( n - 1, 1 ), 1 ) + ... 
    diag( -2 * ones( n - 1, 1 ), - 1 ) + ...
    diag( -1 * ones( n - 2, 1 ), - 2 );
A( n, 1 ) = - 1;
A( 1, n ) = 1;
figure; spy( A )
[ L, U, P ] = lu( A );
figure; spy( L )
figure; spy( U )

n = 1000;
As = sparse( 1 : n, 1 : n, 4 * ones( n, 1 ), n, n ) + ...
     sparse( 1 : n - 1, 2 : n, -1 * ones( n - 1, 1 ), n, n ) + ...
     sparse( 2 : n, 1 : n - 1, -2 * ones( n - 1, 1 ), n, n ) + ...
     sparse( 3 : n, 1 : n - 2, -1 * ones( n - 2, 1 ), n, n );
As( n, 1 ) = - 1;
As( 1, n ) = 1;
A = diag( 4 * ones( n, 1 ), 0 ) + ...
    diag( -1 * ones( n - 1, 1 ), 1 ) + ... 
    diag( -2 * ones( n - 1, 1 ), - 1 ) + ...
    diag( -1 * ones( n - 2, 1 ), - 2 );
A( n, 1 ) = - 1;
A( 1, n ) = 1;
whos A As
%  Name         Size                Bytes  Class     Attributes
%
%  A         1000x1000            8000000  double              
%  As        1000x1000              72104  double    sparse

A = [ 4 -2 -1; -2 7 -4; -1 -4 6 ];
x_ex = ones( 3, 1 ); b = A * x_ex;
R = chol( A );
y = forward_substitutions( R', b );
x = backward_substitutions( R, y );
err = norm( x - x_ex )
% err =
%   9.1551e-16