function x=solve(LR,b)
%*** get size of matrix
n = size(LR,1);

x = zeros(n,1); 
%*** solve Ly=b and save result in b 
for j = 1:n
  for k=1:j-1   
    b(j) = b(j) - LR(j,k) * b(k);
  end
end

%*** solve Rx=y
for j = n:-1:1
  for k=j+1:n 
    b(j) = b(j) - LR(j,k) * x(k);
  end
  x(j) = b(j)/LR(j,j);
end