function [x]=lu_fac(a,b)
% Usage: [x]=lu_fac(a,b)
% Solves the linear system ax = b by lu factorization LUx = b
% Input:
% a = matrix
% b = right hand side vector
% Output:
% x = solution vector

[l,u]=gauss_elim(a); % Calculate LU factorization of a
x = backsub(u,forsub(l,b)); % Solve LUx = b