function [L, U, x] = my_LU_solver(n, a, b)

% a = [a,b];

% -------- LDL' decomposition -----------

L = eye(n,n); U = zeros(n,n); 

 for k = 1:n
	 % Get k-th row of U
    for j=k:n
       U(k,j) = a(k,j);
       for l=1:k-1
	 U(k,j) = U(k,j)-L(k,l)*U(l,j);  
       end
    end

    for i=k+1:n
	  % Get i-th row of L
       L(i,k) = a(i,k);
       for j=1:k-1
	   L(i,k) =  L(i,k) -L(i,j)*U(j,k);
       end
        L(i,k) =  L(i,k)/U(k,k);
    end
 end

 y=zeros(n,1); % Ly = b using FW substitution
 for i=1:n
   y(i) = b(i);
   for j=1:i-1
     y(i) = y(i) - L(i,j)*y(j);
   end
end   

 x= zeros(n,1);  % Ux=y  using BW substitution
 for i=n:-1:1
   x(i) = y(i);
   for j=i+1:n
     x(i) = x(i) - U(i,j)*x(j);
   end
   x(i) = x(i)/U(i,i);
 end