我写了一个Matlab的化学动力学计算程序，可是运行时有问题，那位大侠帮忙指导下，谢谢

function fourlump
clear all
clc
global theta yexp V21 V23 V24 A xx0 M MH x0
k0 = [1.0  5.0  0.5];         % 参数初值
lb = [0  0  0];                   % 参数下限
ub = [+inf  +inf  +inf];    % 参数上限
xx0=[31.423 42.861 6.403 19.315];
M=[93.99229 83.27657 102.8053 116.2428];
NH=[0.2271 0.2271 0.2271 0.0852 0.0852 0.0852 0.1135 0.1135 0.1419 0.1419 0.1419 0.1419 0.1703 0.1703 0.1703];

tspan=[0 1];
ExpData= ...
[1        0.002553        0.002189        0.000629        0.001467
2        0.003354        0.001745        0.000684        0.001973
3        0.002758        0.002128        0.000585        0.001336
4        0.003226        0.003156        0.000672        0.001585
5        0.003495        0.002759        0.000707        0.001658
6        0.003105        0.002974        0.000678        0.001677
7        0.003571        0.001938        0.000791        0.001780
8        0.003073        0.002763        0.000675        0.001594
9        0.002857        0.002761        0.000630        0.001569
10        0.003621        0.001888        0.000679        0.001637
11        0.003228        0.002415        0.000629        0.001522
12        0.003121        0.001553        0.000734        0.001877
13        0.002850        0.002348        0.000636        0.001601
14        0.003406        0.001851        0.000648        0.001570
15        0.003142        0.001809        0.000677        0.001687];
yexp=ExpData(:,2:5);

tspan=[0 1];
V21=[0.8837 1.0897 0.9285 0.8446 0.9119 0.8925 0.9612 0.8992 0.8601 0.9473 0.9115 0.9747 0.8795 0.9365 0.9162]';
V23=[0.8853 0.9914 0.8874 0.8228 0.8501 0.8609 0.9220 0.8669 0.8354 0.8654 0.8592 0.9637 0.8455 0.8638 0.9014]';
V24=[0.7570 0.8668 0.7843 0.7332 0.7551 0.7666 0.7936 0.7678 0.7411 0.7644 0.7638 0.8576 0.7503 0.7658 0.8033]';
theta=[166.8411 250.2616 286.0133 155.4816 248.7706 266.5399 169.4721 593.1523 80.9890 377.9487 566.9230 566.9229 180.9717 361.9435 434.3321];
A=[0.1562 0.1562 0.1562 0.0728 0.0728 0.0728 0.0925 0.0925 0.1106 0.1106 0.1106 0.1106 0.1270 0.1270 0.1270]';
yexp=ExpData(:,2:5);

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunclump,k0,lb,ub,[],tspan,xx0,x0,theta,V21,V23,V24,yexp)
ci = nlparci(k,resid,jacobian)

function f = ObjFunclump(k,tspan,x0,theta,V21,V23,V24,A,yexp)
for i=1:length(theta)
    x0(1)=xx0(1)/((1+2*NH(i))*M(1));
    x0(2)=xx0(2)/((1+2*NH(i))*M(2));
    x0(3)=xx0(3)/((1+2*NH(i))*M(3));
    x0(4)=xx0(4)/((1+2*NH(i))*M(4));
    [t x] = ode45(@FCClump,tspan,x0,[],k,theta(i),V21(i),V23(i),V24(i),A(i));
    y(i,1) = x(end,1);
    y(i,2:4) = x(end,2:4);
end
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f = [f1;f2;f3;f4];

function dxdt=FCClump(t,x,k,theta,V21,V23,V24,A)

dxdt =  ...
[ theta*k(1)*V21*x(2)/(x(1)+x(2)+x(3)+x(4)+A)
-1*theta*(k(1)+k(2)+k(3))*x(2)/(x(1)+x(2)+x(3)+x(4)+A)
  theta*k(2)*V23*x(2)/(x(1)+x(2)+x(3)+x(4)+A)
  theta*k(3)*V24*x(2)/(x(1)+x(2)+x(3)+x(4)+A)];

什么问题