%% subroutine of LU decomposition (solving Tridiagonal System) %%
function [u1,d1,l1] = decomt(u,d,l)
 
n = length(d); % get the length of the diagonal
 
% initialize the output vectors
u1=u;
d1=d;
l1=l;
 
%Perform LU decomposition
d1(1)=d(1);
for i=2:n
    l1(i-1) = l(i-1)/d1(i-1); % Update the lover triangle vector
    d1(i) = d(i) - (l(i-1)/d1(i-1))*u(i-1); %Update the diaginal    
end


