MATLAB调用eval函数时遇到了问题怎么解决

想求解一个非线性方程组，然后设置了一系列符号变量，但是最后符号变量转换求解时报错了，有没有人帮我看看，非常感谢

clear;clc;close all
%%定义方程组
n=3;  %裂缝网络的断裂数目
for i=1:2*n  %定义符号变量
  x(i)=sym(['x',num2str(i)]);
  f(i)=sym(['f',num2str(i)]);
end
%%定义非线性方程
f(1)=x(2);
for i=1:n-1
    f(2*i)=(-1/3*(1/6*((i+1)/n)^3-1/6*(i/n)^3+1/2*((i+1)/n)^2*x(2*i+1)-1/2*(i/n)^2*x(2*i+1)+1/n*x(2*i+2)))^(2/3)*(1/2*(i/n)^2+i/n*x(2*i-1)+x(2*i))-...
        (-1/3*(1/6*(i/n)^3-1/6*((i-1)/n)^3+1/2*(i/n)^2*x(2*i-1)-1/2*((i-1)/n)^2*x(2*i-1)+1/n*x(2*i)))^(2/3)*(1/2*(i/n)^2+i/n*x(2*i+1)+x(2*i+2));
    f(2*i+1)=2*(-1/3*(1/6*((i+1)/n)^3-1/6*(i/n)^3+1/2*((i+1)/n)^2*x(2*i+1)-1/2*(i/n)^2*x(2*i+1)+1/n*x(2*i+2)))^(1/3)*(i/n+x(2*i+1))-...
        (-1/3*(1/6*(i/n)^3-1/6*((i-1)/n)^3+1/2*(i/n)^2*x(2*i-1)-1/2*((i-1)/n)^2*x(2*i-1)+1/n*x(2*i)))^(1/3)*(i/n+x(2*i-1));
end
f(2*n)=x(2*n)+1;
f=[f(1);f(2);f(3);f(4);f(5);f(6)];  %方程组成方程组
%%计算雅克比矩阵
J=[diff(f(1),x(1)) diff(f(1),x(2)) diff(f(1),x(3)) diff(f(1),x(4)) diff(f(1),x(5)) diff(f(1),x(6))
    diff(f(2),x(1)) diff(f(2),x(2)) diff(f(2),x(3)) diff(f(2),x(4)) diff(f(2),x(5)) diff(f(2),x(6))
    diff(f(3),x(1)) diff(f(3),x(2)) diff(f(3),x(3)) diff(f(3),x(4)) diff(f(3),x(5)) diff(f(3),x(6))
    diff(f(4),x(1)) diff(f(4),x(2)) diff(f(4),x(3)) diff(f(4),x(4)) diff(f(4),x(5)) diff(f(4),x(6))
    diff(f(5),x(1)) diff(f(5),x(2)) diff(f(5),x(3)) diff(f(5),x(4)) diff(f(5),x(5)) diff(f(5),x(6))
    diff(f(6),x(1)) diff(f(6),x(2)) diff(f(6),x(3)) diff(f(6),x(4)) diff(f(6),x(5)) diff(f(6),x(6))];
%%定义牛顿迭代法求方程组的解
n=1;  %迭代次数
y0=[-5;-5;-5;-5;-5;-5];  %迭代初值
E=1; 
while E > 1e-4  %如果精度不满足要求就一直迭代
    x(1)=y0(1); %给变量赋值
    x(2)=y0(2);
    x(3)=y0(3);
    x(4)=y0(4);
    x(5)=y0(5);
    x(6)=y0(6);
    y1=y0-inv(eval(J))*eval(f); %迭代计算，用x1替代x0
    E=max(abs((y1-y0)./y1));    %计算求解的精度
    y0=y1;    %将当前的解x1赋值给x0，为下一次迭代做准备
    n=n+1;   %迭代次数加一
end
%%显示方程组的解，并将求解结果带回方程组，验算求解结果
x(1)=y1(1)
x(2)=y1(2)
x(3)=y1(3)
x(4)=y1(4)
x(5)=y1(5)
x(6)=y1(6)
eval(f)

img

eval的形参是一个表达式，你不能传一个矩阵进去吧？

img