Numerical Methods HW 2

function A=lugauss(A)
%LUGAUSS LU factorization without pivoting
% A = LUGAUSS(A) stores an upper triangular matrix in
% the upper triangular part of A and a lower triangular
% matrix in the strictly lower part of A (the diagonal
% elements of L are 1).
[n,m]=size(A);
if n ~= m; error('A is not a square matrix'); else
 for k = 1:n-1
  for i = k+1:n
   A(i,k) = A(i,k)/A(k,k);
   if A(k,k) == 0, error('Null diagonal element'); end
   j = [k+1:n]; A(i,j) = A(i,j) - A(i,k)*A(k,j);
  end
 end
end

for n=1:15
    % Create A
    A = eye(n);
    for i = 1:n
        for j = 1:n
            A(i,j) = 1/(i+j-1);
        end
    end

    % Create b
    b = transpose(zeros(1,n));
    b(1) = 1;

    % LU-decomposition
    for k = 1:n-1
        for i = k+1:n
            A(i,k) = A(i,k)/A(k,k);
            if A(k,k) == 0, error('Null diagonal element'); end
            j = [k+1:n]; A(i,j) = A(i,j) - A(i,k)*A(k,j);
        end
    end
    L = eye(n) + tril(A,-1);
    U = triu(A);

    % Create y
    y = transpose(zeros(1,n));

    % Forward substitution
    y(1) = b(1);
    for i = 2:n; y = [y; b(i)-A(i,1:i-1)*y(1:i-1)]; end

    % Backward substitution
    x(n)=y(n)/A(n,n);
    for i = 3:-1:1
        x(i) = (y(i)-A(i,i+1:n)*x(i+1:n)')/A(i,i);
    end

    A*x - b

    fprintf('\nn=%i\n',n) % Print data we care about.
    K=cond(A)
    alpha=K*(norm(b-A*x))/(norm(b))
end