function [L,U,alpha,gamma]=LU_decomposition(a,b,c,JL)

%a is the main diagonal 
%b is the lower diagonal
%c is the upper diagonal

%%Start by obtaining the first row
alpha=zeros(JL,1);      gamma=zeros(JL-1,1);

alpha(1)=a(1);      
gamma(1)=c(1)/alpha(1);

for i=2: JL-1
    alpha(i)=a(i)-b(i-1)*gamma(i-1);
    gamma(i)=c(i)/alpha(i);
    
end

alpha(JL)=a(JL)-b(JL-1)*gamma(JL-1);

L=diag(alpha)+diag(b,-1);
U=diag(ones(JL,1))+diag(gamma,1);