MATLAB: Newton-Raphson Method with Jacobian

jacobiannewtonnewton raphsonraphson

I have a problem with this program, a finite value vector is not returned despite the system having a solution.
Using function fsolve the result is Xeq3 = [0.6875 0.6346 0.9411], while using the function my_newton2 Xeq3 = [NaN NaN NaN].
I think the problem is in the function declaration but I was unable to find a solution. Please help me to find the error in the code.
function es()
    KI = 1;
    k1 = 1.2;
    k2 = 1.1;
    k3 = 1.3;
    tf = 10;
    X0 = [2, 1.25, 0.75];
%==================================
    format long
    syms x y z;
    f(x,y,z) = [(k3 / (z + y) - k1 * x), (k1 * x - k2 * 0.75), (k2 * 0.75 - k3 * y)];
    j(x,y,z) = jacobian(f, [x, y, z]);
    j1(x,y,z) = inv(j);
    
    Xeq3 = my_newton2(f, j1, [X0(1), X0(3), KI])
end
%===================================
function xval = my_newton2(f, j, x0)
    tol = 1e-9; %tollerance
    while true 
        % Calculate the next xval value
        h = -f(x0(1),x0(2),x0(3))*j(x0(1),x0(2),x0(3));
        xval = double(x0) + double(h);
        
	%Check the result
        exit = true;
        for k=1:length(x0)
            if abs(x0(k)-xval(k))>tol
                exit=false;
             end
        end
        if exit
            break;
        end
        x0 = xval;
    end
end

Best Answer

In Newton loops you must evaluate your f and j and currently guessed iterate, so your line
h = -f(x0(1),x0(2),x0(3))*j(x0(1),x0(2),x0(3));
should rather look like (using your notation)
h = -f(xval(1),xval(2),xval(3))*j(xval(1),xval(2),xval(3));
so you will need a statement at the top to initialize xval to x0
Related Question