剛入門matlab的小白，最近想做反應的動力學，按照B站up主的視頻自己寫了一段代碼，但是運行總是出問題：

錯誤使用 odearguments (第 93 行)；FUNC 必須返回列向量，出錯 ode45 (第 115 行)，odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
出錯 Kinetics>fun (第 67 行)，[t,x]=ode45(@func,tspan,x0,[],k)；
出錯 lsqnonlin (第 218 行)，initVals.F = feval(funfcn{3},xCurrent,varargin{:});出錯 Kinetics (第 18 行)，lsqnonlin(@fun,k0,lb,ub,[],yexp);%非線性最小二乘法。原因:Failure in initial objective function evaluation. LSQNONLIN cannot continue.

下面是我寫的代碼，讀取的Excel表格里有5列*7行的實驗數(shù)據(jù)，勞煩大佬幫我瞅瞅哪里需要改動，萬分感謝。

function Kinetics
%反應一：A+B=C+M
%r=k*XA*XB-K*XC*XM
%反應二：A+C=D+M
%r=K*XA*XC-K*XD*XM
%反應三：B+C=E+M
%r=K*XB*XC-K*XE*XM
%XM=0.175
clc
clear all;
global a b
tspan=[0.5 1 4 6 8 12 16];
yexp=xlsread('reaction.xls');
k0=[0.1 0.01 0.01 0.001 0.001 0.001];%參數(shù)初值
lb=[0 0 0 0 0 0];%下邊界
ub=[+inf +inf +inf +inf +inf +inf];%上邊界
[k,resnorm,residual,exitflag,output,lambda,jacobian]=...
    lsqnonlin(@fun,k0,lb,ub,[],yexp);%非線性最小二乘法
tspan=[0.5 1 4 6 8 12 16];
a=1;
b=a+6;
x0=yexp(a,;%積分初值
[t,x]=ode45(@func,tspan,x0,[],k);
t1=linspace(0.5,16,200);
ya1=spline(t,x(:,1),t1);%動力學計算得到的點進行樣條插值
ya2=spline(t,x(:,2),t1);
ya3=spline(t,x(:,3),t1);
ya4=spline(t,x(:,4),t1);
ya5=spline(t,x(:,5),t1);
for m=1:7
    for n=1:5
        yy(a+m-1,n)=x(m,n);%每一次的值存入yy矩陣
    end
end
figure(1)
plot(tspan,yexp(a:b,1),'k^',t1,ya1,'k-',tspan,yexp(a:b,2),'ro',t1,ya2,'r-',tspan,yexp(a:b,3),'bd',t1,ya3,'b-',...
tspan,yexp(a:b,4),'g*',t1,ya4,'g-',tspan,yexp(a:b,5),'yp',t1,ya5,'y-');
legend('','A濃度','','B濃度','','C濃度','','D濃度','','E濃度');
xlabel('t(h)');ylabel('濃度(mol/L)');title('170℃ 0.1wt%催化劑');
t1=linspace(0.5,16,200);
z1=spline(t,yy(1:7,1),t1);
h1=spline(t,yy(1:7,2),t1);
s1=spline(t,yy(1:7,3),t1);
b1=spline(t,yy(1:7,4),t1);
u1=spline(t,yy(1:7,5),t1);
xlswrite('result.xls',[t1' z1' h1' s1' b1' u1'],'sheet1');
xlswrite('result.xls',residual,'sheet2');
Ne = length(yexp(:,2));     %模型適定性判別
Np = length(k);
[rho2,F] = rho2_F(k,yexp,resnorm,Ne,Np);
ci=nlparci(k,residual,jacobian)
fprintf('\t k1,0=%.1f ± %.4f\n',k(1),ci(1,2)-k(1));
fprintf('\t k2,0=%.1f ± %.4f\n',k(2),ci(2,2)-k(2));
fprintf('\t k3,0=%.1f ± %.4f\n',k(3),ci(3,2)-k(3));
fprintf('\t k4,0=%.1f ± %.4f\n',k(4),ci(4,2)-k(4));
fprintf('\t k5,0=%.1f ± %.4f\n',k(5),ci(5,2)-k(5));
fprintf('\t 殘差平方和：%.3f\n',resnorm)
fprintf('\t 實驗點數(shù)和自由度分別為 Ne = %d和 Np = %d\n',Ne,Np)
fprintf('\t 決定性指標ρ^2: %.4f\n',rho2)
fprintf('\t F比: %.3f\n\n',F)
%=================================================================================
function f=fun(k,yexp)
f=[];
tspan=[0.5 1 4 6 8 12 16];
a=1;
x0=yexp(a,;
[t,x]=ode45(@func,tspan,x0,[],k)
d=a+6;
yc1=x(:,1);
yc2=x(:,2);
yc3=x(:,3);
yc4=x(:,4);
yc5=x(:,5);
f11=yexp(a:d,1)-yc1;
f12=yexp(a:d,2)-yc2;
f13=yexp(a:d,3)-yc3;
f14=yexp(a:d,4)-yc4;
f15=yexp(a:d,5)-yc5;
ff=[f11 f12 f13 f14 f15];
f=[f;ff];
%=================================================================================
function dxdt=func(t,x,k)
r1=-k(1)*x(1)*x(2)-k(2)*x(1)*x(3)+k(4)*x(3)*0.175+k(5)*x(4)*0.175;
r2=-k(1)*x(1)*x(2)-k(3)*x(2)*x(3)+k(4)*x(3)*0.175+k(6)*x(5)*0.175;
r3=k(1)*x(1)*x(2)+k(5)*x(4)*0.175+k(6)*x(5)*0.175-k(2)*x(1)*x(3)-k(3)*x(2)*x(3)-k(4)*x(3)*0.175;
r4=k(2)*x(1)*x(3)-k(5)*x(4)*0.175;
r5=k(3)*x(2)*x(3)-k(6)*x(5)*0.175;
dxdt=[r1 r2 r3 r4 r5]
%=================================================================================
function [rho2,F] = rho2_F(k,yexp,s,Ne,Np)
y=yexp.^2;
sy = sum(y();
rho2 = 1 - s/sy;              %rho2: 決定性指標
F = (sy - s)*(Ne-Np)/(Np*s);  %F：F比 返回小木蟲查看更多