方程组为：x1^3+x1^2-x2+x3-13=0; x2^3+x2^2-x1+2x3-20=0; x3^3+x3^2+2x1-2x2-18=0，其精确解是（2,2,3）参考初值为（1，1，2）
注：x1,x2,x3为未知数，最高三次，而且是方程组。
function [r,n]=mulNewton(F,x0,eps)
if nargin==2
eps=1.0e-4;
end
x0 = transpose(x0);
Fx = subs(F,findsym(F),x0);
var = findsym(F);
dF = jacobian(F);
dFx = subs(dF,findsym(dF),x0);
r=x0-inv(dFx)*Fx;
n=1;
tol=1;
while tol>eps
x0=r;
Fx = subs(F,findsym(F),x0);
dFx = subs(dF,findsym(dF),x0);
r=x0-inv(dFx)*Fx; %核心迭代公式
tol=norm(r-x0);
n=n+1;
if(n>100000) %迭代步数控制
disp('迭代步数太多，可能不收敛！');
return;
end
end

供参考