function [ y ] = gausselim( A, b )
%UNTITLED Solves the system Ay=b for y
%   Input is n x n matrix A and n x 1 vector b

n=length(A);
q=1:n;
if det(A)==0
    error('det(A)=0')
else
    for k=1:n-1
%Find pivot point and keep track of what row and column in came from
[colmax, rowindices]=max(abs(A(k:n,k:n)));
[biggest, colindex]=max(colmax);
row=rowindices(colindex)+k-1;
col=colindex+k-1;

%Switch the rows and columns
A([k, row],:) = A([row, k],:);
b([k, row],:) = b([row, k],:);
A(:,[k, col]) = A(:,[col, k]);
q([k, col]) = q([col, k]);

%Gauss elimination
for i = 1:n-1
  for j = i+1:n
    m = A(j,i)/A(i,i);
    A(j,:) = A(j,:) - m*A(i,:);
    b(j) = b(j) - m*b(i);
  end
end

    end
%back substitution
x=zeros(1,n);
x(n)=b(n)/A(n,n);
for l=1:n-1 
    x(n-l)=(b(n-l)-(dot(A(n-l,:),x)))/A(n-l,n-l);
end

%back to right order
y=zeros(n,1);
for m=1:n
    for p=1:n
        if q(p)==m
            y(m)=x(p);
        end
    end
end

end
end   


