求助，这是求解非线性方程组的代码，改成解有多个变量的，何解？

基于我现在手头的学习任务，我想把这个代码改造成可以求解含有多个乃至于成百上千个变量x1,...,xn的方程组的代码，方程组为Cx-f=0，其中C为含有变量x1,...,xn的矩阵，f为含有变量x1,...,xn的向量，但问题是，多个变量定义起来很麻烦，如何实现“syms x1 x2 ... xn;x=[x1,x2,...,xn];”？
以下是使用牛顿迭代法求解二元非线性方程组的例题代码（初值[x y]=[0 0 ]，精度为1e-5）:
①
function [ f,x ] = fun
syms x1 x2;
x=[x1,x2];%x为函数变量
f1=x1^2-10*x1+x2^2+8;
f2=x1*x2^2+x1-10*x2+8;
f=[f1;f2];%f为所定义的函数
end
②
function[X,time]=newton(x0,eps,N)
%x0为初始值
%eps为允许误差
%N为最大迭代次数
[f,x]=fun;
n=length(x);
for i=1:n
df(:,i)=diff(f,x(i));%按列依次求偏导
end
time=0;
while(time<=N)%进行迭代过程
fx=f;
dfx=df;
for i=1:n
fx=subs(fx,x(i),x0(i));
end
for i=1:n
dfx=subs(dfx,x(i),x0(i));%将初始值赋入，计算函数值
end
fx=eval(fx);
dfx=eval(dfx);
X=x0-dfx\fx;
time=time+1;
error=norm(X-x0);%计算误差
if error<eps %与允许误差进行比较
break;
end
x0=X;
end
if time>N
warning('超出迭代的次数');
end
③
[X,time]=newton([0;0],1e-5,100)