[MATLAB数学相关] 问题:matlab实现牛顿迭代法求解非线性方程

[matlab0407]

最近在研究weibull分布的参数估计，其中有一个环节为matlab编程利用最大似然函数求解估计参数，在此就不可避免的要面临求解非线性方程。
本人在网上搜索多个文献，其中被大家引用最多的一个例子就是"matlab实现牛顿迭代法求解非线性方程"。通过试着运行作者的matlab code，产生了如下一些疑问，请各位老师帮忙解答，不胜感谢。

clear
clc

syms x1 x2 x3 % 定义三个未知数

f1=3*x1-cos(x2*x3)-1/2;
f2=x1^2-81*(x2+0.1)^2+sin(x3)+1.06;
f3=exp(-x1*x2)+20*x3+(10*pi-3)/3; % 三个方程组成一个方程组并且f1=0;f2=0;f3=0

f=[f1 f2 f3];
df=[diff(f,x1);diff(f,x2);diff(f,x3)]; % 分别对x1 x2 x3 三个未知数求偏导

x0=[0.1 0.1 -0.1];
eps=0.00001;
N=20; % 设置初始值 x1 x2 x3 的初始值分别为0.1 0.1 -0.1，精确度0.00001 迭代运算次数设置为最多20次

% 求解
for i=1:N; % 运算次数

f=subs(f,{'x1' 'x2' 'x3'},{x0(1) x0(2) x0(3)}); % 为subs函数的赋值运算 1x3的矩阵
df=subs(df,{'x1' 'x2' 'x3'},{x0(1) x0(2) x0(3)});% 为subs函数的赋值运算 3x3的矩阵

x=x0-f/df; % 得到的一组新数值 1x3的矩阵 理解为 新的x0=[新x0(1) 新x0(2) 新x0(3)]

if norm(x-x0)<eps
break;
end % 条件判断 得出的新值是否满足精确度的要求 可以求x矩阵元素平方和的开根

x0=x;% 当条件满足 输出求得的x0

end

运行上述程式出现的问题如下
当第一次循环初始值为x0=[0.1 0.1 -0.1]代入，得到第一次循环的计算结果 f/df得到一个结果 进而x得到一个结果
但是当第二次循环 f/df的值就已经保持不变 仅有x和x0在循环计算
当然这不是我们想要
请各位老师指点
PS 有没有更好的迭代的在matlab里的实现方法可以更好的缩短结算时间？