Newton-Raphson Method with Jacobian

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

Answer 1

[J. Alex Lee]

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

[Morena Fabozzi]

I update x0 at the end of each loop. However, I corrected what you told me but I still have the same result.

[J. Alex Lee]

By the way I didn't see that at the bottom of your newton method you overwrote x0 with xval, which achieves the same effect as what I suggested; so clearly that wasn't your problem.

Did you actually try with fsolve to get the solution Xeq3 = [0.6875 0.6346 0.9411]? What happens if you start your newton method with this "correct" solution?

Can you monitor the value of h from the very first iteration?

[Morena Fabozzi]

Yes I tried to use fsolve for the corresponding function handler and the result is correct. Using the method with x0 = [0.6875 0.6346 0.9411] the result is correct. I can monitor h but I don't know what values ​​it must have. What do you think the problem is? Newton's method or functions f and j?

[J. Alex Lee]

If your newton's method works fine when your initial guess is close to the actual solution, then your newton's method is taking steps out of the region of convergence, and depending on your equations may be stepping into domains where your function values are NaN. This is just a limitation of newton's method, which requires reasonably good initial guesses. fsolve may be more robust to this issue.

[Morena Fabozzi]

It works. Thanks