clear all
 n=input('n=')
 onOnes = ones(n, 1) ;
 % a = diag(2 * nOnes, 0) - diag(nOnes(1:n-1), -1) - diag(nOnes(1:n-1), 1)

 % Crout LU decomposotion; Demonstration only. We usually  do not form the matrix

 d = 1:n; d=100*d;  a=1:n-1; b=-2*a; 
 A = diag(d,0)+ diag(a,-1) + diag(b,1);

 dn = d; an =a;
 for i=2:n
   an(i-1) = an(i-1)/dn(i-1);
   dn(i) = dn(i) - an(i-1)*b(i-1);
 end

xe = ones(n,1); b1=A*xe;

[L, U, x] = my_LU_solver(n, A, b1)


 
% Can you solve A x = b1 using only an, dn, and b?