请问MATLAB中利用牛拉法解含参数的方程如何求解呢？

目的是输出X有一定范围，根据这个范围反推出f(x)系数的取值范围。我的思路是将变系数作为字母a传入函数，最后解得的X结果中含有字母a，再根据限制范围计算字母a的范围。该常数只存在于f()的个别行中
尝试使用sysm 但是报错了，其中求解f(x)使用的函数:
disp(norm(f(x(1,:))));
for n=2:length(x)
x(n,:)=x(n-1,:)-((J(x(n-1,:))^(-1))*f(x(n-1,:)))';
disp(norm(f(x(n,:))));
end
感谢指导！
报错显示：
从 sym 转换为 double 时出现以下错误:
Unable to convert expression into double array.

出错 try_d (line 68)
x(n,:)=x(n-1,:)-((Jp(x(n-1,:))^(-1))*fp(x(n-1,:)))';

你好同学，建议用牛顿拉夫逊方法时不要使用符号，那样太拖慢速度了，而且符号变量与双精度实型变量之间转换特别麻烦。给你举个简单的例子你看看
如果用符号，尝试一下：

x(n,:)=x(n-1,:)-((J(x(n-1,:))^(-1))*f(x(n-1,:)))';

你可以参考我之前的N-R系列


clear
clc

V1=1.01;

ZL12=0.02+i0.01;
ZL23=0.03+i0.01;
Zload2=9.8+i1.96;
Sload3=0.2+i0.05;

YL12=1/ZL12;
YL23=1/ZL23;
Yload2=1/Zload2;
Pload3=real(Sload3);
Qload3=imag(Sload3);
sysm Qp=0.2;
%这里求导后的J里面不含该参数，但是方程第一二行是含有的。所以我在想迭代的时候会不会出错
%另外，想要根据迭代的最后结果V在0.8-1.2之间来求得Pp的取值范围，不晓得我的问题表达清楚了没有
%sysm Sp=0.8;
%Pp=real(Sp);
%Qp=imag(Sp);

Y=[YL12,-YL12,0;
-YL12,YL12+YL23+Yload2,-YL23;
0,-YL23,YL23];

G12=real(Y(1,2));
B12=imag(Y(1,2));

G22=real(Y(2,2));
B22=imag(Y(2,2));

G23=real(Y(2,3));
B23=imag(Y(2,3));

G33=real(Y(3,3));
B33=imag(Y(3,3));

f=@(V2,V3,theta2,theta3) [V1*V2*(G12*cos(theta2)+B12*sin(theta2))+V2*V2*G22+V2*V3*(G23*cos(theta2-theta3)+B23*sin(theta2-theta3))+Pp; V2*V3*(G23*cos(theta3-theta2)+B23*sin(theta3-theta2))+V3*V3*G33+Pload3-Pp; V1*V2*(G12*sin(theta2)-B12*cos(theta2))-V2*V2*B22+V2*V3*(G23*sin(theta2-theta3)-B23*cos(theta2-theta3)); V2*V3(G23*sin(theta3-theta2)-B23*cos(theta3-theta2))-V3*V3*B33+Qload3];

J=@(V2,V3,theta2,theta3) [V1*(G12*cos(theta2)+B12*sin(theta2))+2*V2*G22+V3*(G23*cos(theta2-theta3)+B23*sin(theta2-theta3)),V2*(G23*cos(theta2-theta3)+B23*sin(theta2-theta3)),V1*V2*(-G12*sin(theta2)+B12*cos(theta2))+V2*(-G23*sin(theta2-theta3)+B23*cos(theta2-theta3)),V2*V3*(G23*sin(theta2-theta3)-B23*cos(theta2-theta3));V3*(G23*cos(theta3-theta2)+B23*sin(theta3-theta2)),V2*(G23*cos(theta3-theta2)+B23*sin(theta3-theta2))+2*V3*G33,V2*V3*(G23*sin(theta3-theta2)-B23*cos(theta3-theta2)),V2*V3*(-G23*sin(theta3-theta2)+B23*cos(theta3-theta2));V2*V1*(G12*sin(theta2)-B12*cos(theta2))-2*V2*B22+V3*(G23*sin(theta2-theta3)-B23*cos(theta2-theta3)),V2*(G23*sin(theta2-theta3)-B23*cos(theta2-theta3)),V1*V2*(G12*cos(theta2)+B12*sin(theta2))+V2*V3*(G23*cos(theta2-theta3)+B23*sin(theta2-theta3)),V2*V3*(-G23*cos(theta2-theta3)-B23*sin(theta2-theta3));V3*(G23*sin(theta3-theta2)-B23*cos(theta3-theta2)),V2*(G23*sin(theta3-theta2)-B23*cos(theta3-theta2))-2*V3*B33,V2*V3*(-G23*cos(theta3-theta2)-B23*sin(theta3-theta2)),V2*V3*(G23*cos(theta3-theta2)+B23*sin(theta3-theta2))];

fp=@(x) f(x(1),x(2),x(3),x(4));
Jp=@(x) J(x(1),x(2),x(3),x(4));

x=zeros(3,4);
x(1,:)=[1.01,1.01,0,0];
disp(norm(fp(x(1,:))));
for n=2:length(x)
x(n,:)=x(n-1,:)-((Jp(x(n-1,:))^(-1))*fp(x(n-1,:)))';
disp(norm(fp(x(n,:))));
end
disp('Solution:');
disp(x(n,:));