matlab牛顿迭代法解非线性方程求解

function f=fun(x);
%定义非线性方程组如下
%变量x1 x2
%函数f1 f2
syms x1 x2
f1 = sqrt((x1-4)^2 + x2^2)-sqrt(17);
f2 = sqrt(x1^2 + (x2-4)^2)-5;
f=[f1 f2];

function df=dfun(x);
f=fun(x);
df=[diff(f,'x1');diff(f,'x2')]; %雅克比矩阵

 clear;clc
 format;
 x0=[0 0];   % 迭代初始值
 eps = 0.00001;  % 定位精度要求
 for i = 1:10
     f = double(subs(fun(x0),{'x1' 'x2'},{x0(1) x0(2)}));
     df = double(subs(dfun(x0),{'x1' 'x2'},{x0(1) x0(2)}));  % 得到雅克比矩阵
     x = x0 - f/df;
     if(abs(x-x0) < eps)
         break;
     end
     x0 = x; % 更新迭代结果
 end
 disp('定位坐标：');
 x
 disp('迭代次数：');
 i