What is the Matlab code for Newton Raphson method to solve the polynomial system? I  want to approximate the zeros of the following system of polynomial equations using the Newton-Raphson method:
\begin{align}
f(x,y) &= x + \frac13 y^9 + \frac19 x^{243} + \frac{1}{27} y^{2187} =0\\ g(x,y) &= y + \frac13 x^{27} + \frac19 y^{243} + \frac{1}{27} x^{6561} =0
\end{align}
I have used the following Matlab code with initial guess $(-1,-1)$ and with error tolerance $0.01$:
close all; clc,clear
% Newton Raphson solution of two nonlinear algebraic equations
xy = [-1 -1]; % initial guesses
iter=0;
maxiter=100;
xy_N = 0.1;
TOL = 0.01;
error1 = [1, 1];
f = zeros(iter, 2);
error = zeros(iter, 2);
% begin iteration
while xy_N>TOL
iter=iter+1;
x = xy(1);
y = xy(2);
% calculate the functions
f(iter,:) = [x+(1/3)*y^9+(1/9)*x^243+(1/27)*y^2187;
y+(1/3)*x^27+(1/9)*y^243+(1/27)*x^6561];
% calculate the Jacobian
J= [1+27*x^242,  3*y^8+81*y^2186; 9*x^26+24*x^6560, 1+27*y^242];
xy_new = xy - ((J)\f(iter,:)')';
error1=abs(xy_new-xy);
error(iter,:)=error1;
% Calculate norm:
xy_N = norm(f(iter,:));
XYN(iter) = xy_N;
xy=xy_new;
%iter=iter+1;
if (iter > maxiter)
fprintf('Did not converge! \n', maxiter);
break
end
end
if iter<maxiter
f(iter+1:end,:)=[];
end
% Print results
fprintf('Number of iterations is:  %f \n ', iter)
fprintf('Final values of [x, y] =   %f  %f \n', xy(end, :));
fprintf('Final values of EQN:  %f  %f \n', f(end, :));
fprintf('Final values of ERROR:  %f  %f \n', error(end, :));

But it is giving error. The following display:
Warning: Matrix is singular to working precision. Warning: Matrix is singular, close to singular or badly scaled.
Results may be inaccurate. RCOND = NaN.
Number of iterations is:  4.000000
Final values of [x, y] =   NaN  NaN
Final values of EQN:  NaN  NaN
Final values of ERROR:  NaN  NaN

Could you please help me fix the code to ensure that convergence is attained? You can suggest more flexible code to solve my system.

Edit:
Since the Jacobian is not invertible, I am adding one more term to each function in order to check my program as follows:
\begin{align}
f(x,y) &= x + \frac13 y^9 + \frac19 x^{243} + \frac{1}{27} y^{2187} +\frac{1}{81}x^{59049}=0\\ g(x,y) &= y + \frac13 x^{27} + \frac19 y^{243} + \frac{1}{27} x^{6561} +\frac{1}{81}y^{59049}=0
\end{align}
I hope now the Jacobian will be non-singular.
 A: Problem seems to be that the Jacobian is not invertible in the fourth iteration. This is because I think the NR-method wants to converge to the solution $x:= (-1.0005, 1.0014)^\top$ but $J(x)$ is not invertible.
The convergence results for the NR-method usually all apply to non-degenerate solutions $y$, i.e. $\det{J(y)} \neq 0$ (of course you cannot check that in advance).
One suggestion that comes to mind is just to store the iterations in which $\texttt{J}$ is invertible and use that old $\texttt{J}$ when it is not invertible in the current iteration (so called frozen NR-method). Even though I do not know for certain if this converges, chance is certainly higher.
(Fun fact: This would actually converge locally - though not as fast as normal NR-method - if $J(x)$ was invertible)
Another one is using the pseudo-inverse of $J$ instead of the $\texttt{\\}$ operator.
A: I have simulated your code step-by-step. On the third iteration, it blowed up.
>>test_math
Number of iterations is:  1.000000 
 Final values of [x, y] =   -0.958268  -0.996274 
Final values of EQN:  -1.481481  -1.481481 
Final values of ERROR:  0.041732  0.003726 
Number of iterations is:  2.000000 
 Final values of [x, y] =   2.610643  -1.777055 
Final values of EQN:  -1.280603  -1.280603 
Final values of ERROR:  3.568911  0.780782 
Warning: Matrix is singular to working precision. 
 In test_math (line 21) 
Number of iterations is:  3.000000 
 Final values of [x, y] =   NaN  NaN 
Final values of EQN:  -Inf  -Inf 
Final values of ERROR:  NaN  NaN 

I have investigated the result and I found out that on the 3rd iteration, the Jacobian matrix has inf value meaning there can be division by zero.
1st iteration
J =

    28    84
    33    28

2nd iteration
J =

    1.0009    2.9349
    2.9710   11.9402

3rd iteration
J =
1.0e+102*1.9243       Inf
   Inf    0.0000

This is probably happen when $x$ or $y$ approaches zero and with very huge exponent.
