clear()
tol = 10e-10
MAX_ITERATIONS = 1000

function [R]=newtonEigen(A)
    R = []
    err = []
    iter = 1
    if size(A,1) ~= size(A,2) then
        printf(" ERROR: Input matrix must be square!\n")
        return
    end
    x = int(rand(size(A,1),1)*10)
    x = x/norm(x)
    lambda = x'*A*x
    n = length(A)
    while iter ~= MAX_ITERATIONS
        //printf("lambda = %f\n",lambda)
        J = [[A -x]; 2*x' 0]
        for i = 1:3
            J(i,i) = J(i,i) - lambda
        end
        B = [A*x-lambda*x ; x'*x-1]
        s = GaussianElim(J,-B)
        x = x+s(1:size(A,1))
        lambda = lambda + s(size(A,1)+1)

        err = [err,norm(s,'inf')]
        if (err($) < tol) then
            q = log10((err($)/err($-1))) / log10((err($-1)/err($-2)))
            printf(" iterations = %d\n q(order of convergence) = %f\n",iter,q)
            disp(lambda,"lambda = ")
            R = x
            return
        end
        iter = iter + 1
    end
endfunction

function [G] = GaussianElim(A, b)
    G = []
    if det(A) == 0 then //check the determinant first
        error("Error: The matrix is linearly dependent")
        return
    elseif size(A,1) ~= size(b,1) //check for size consistency
        error("Error: Check your matrix sizes")
        return
    end
    M = [A, b] //augment the matrices
    rowSize = size(M,1)
    colSize = size(M,2)
    for i = 1:1:rowSize
        //find the pivot
        [pivotValue, pivotIndex] = max(abs(M(i:rowSize,i)))
        pivotIndex = pivotIndex + (i - 1)
        //swap
        tempRow = M(i,:)
        M(i,:) = M(pivotIndex,:)
        M(pivotIndex,:) = tempRow
        //perform ero2 first
        M(i,:) = M(i, :) / M(i,i)
        //i-zero mo na
        for j = i+1:1:rowSize
            mul = M(j,i)
            M(j,:) = M(j,:) - (M(i, :) * mul)
        end
     end
     //start back substitution
         for i = rowSize - 1:-1:1
             for j = i + 1:rowSize
                 M(i,$) = M(i,$) - (M(i,j) * M(j,$))
             end
             M(i,$) = M(i,$) / M(i,i)
         end
     G = M(1:rowSize,$)
endfunction

A = [4,-0.5,0; 0.6,5,-0.6; 0, 0.5 3]
B = [2.9766 0.3645 0.4198 1.1159;
    0.3945 2.7328 -0.3097 0.1123;
    0.4198 -0.3097 2.5675 0.6079;
    1.1159 0.1129 0.6079 1.7231]
disp(A, "A =")
disp(newtonEigen(B),"Eigenvalue of A using newtons =")