function X = gaussj(A,B)
% X = gauss(A,B)
% A is an nxn non-singular matrix, input.
% B is a nx1 vector, input.
% X is an nx1 vector, output.
% Trivial pivoting is used.
[n n] = size(A);
A = [A';B']';
X = zeros(n,1);
for p = 1:n,
for k = [1:p-1,p+1:n],
if A(p,p)==0, break, end
mult = A(k,p)/A(p,p);
A(k,:) = A(k,:) - mult*A(p,:);
end
end
X = A(:,n+1)./diag(A);