function [PHI]=LU_solve(alpha,b,gamma,f,JL)

%%First solve for PHI*==PHI_star where L*PHI_Star=f--------FORWARD SWEEP
PHI=zeros(JL,1);
PHI_STAR=zeros(JL,1);

PHI_STAR(1)=f(1)/alpha(1);

for i=2: JL
    PHI_STAR(i)=(f(i)-b(i-1)*PHI_STAR(i-1))/alpha(i);
end

%%Now solve for PHI since U*PHI_STAR=PHI-----------------BACKWARD SWEEP


PHI(JL)=PHI_STAR(JL);

for j=JL-1:-1:1
    PHI(j)=PHI_STAR(j)-gamma(j)*PHI(j+1);
end