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

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

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

for i=0: JL