function [N,x] = fixpoint (x0)
% Function fixpoint 
% In: Start value
% Out: number of iteration steps, vector of solution
N=1;
xCur = x0;
xOld = x0;
epsilon = 10^(-8);
A = [3,4;1,6];

function [y] = f(x)
    yTemp = (A - x(3)*eye(2))*[x(1);x(2)];
    lambda = x(1)^2+x(2)^2-1;
    y  = [ yTemp ; lambda];
end

function [y] = jacobi(x)
    y1 = [A(1,1)-x(3) , A(1,2)      , -x(1)];
    y2 = [A(2,1)      , A(2,2)-x(3) , -x(2)];
    y3 = [2*x(1)      , 2*x(2)      , 0];
    y  = [y1;y2;y3];
end

function [y] = fixpoint_next(x)
    y = x - inv(jacobi(x0)) * f(x);
end

xCur = fixpoint_next(xOld);
x(:,N) = xCur;
while(norm(xCur-xOld) > epsilon && N < 60)
    N = N + 1;
    xOld = xCur;
    xCur = fixpoint_next(xOld);
    x(:,N) = xCur;
end
end