求助，很着急,求matlab高手帮帮忙吧

我想解一个微分方程组，主程序是这样的
a=0.000000012;
b=0.000000024;
M=3000;
h=(b-a)/M;
t=zeros(1,M+1);
x=zeros(1,M+1);
y=zeros(1,M+1);
z=zeros(1,M+1);
t=a:h:b;
x(1)=6.4164e+010;
y(1)=4372.7;
z(1)=2.2129e+024;
for j=1:M

kx1=h*feval('fx',t(j),x(j),y(j),z(j));
ky1=h*feval('fy',t(j),x(j),y(j),z(j));
kz1=h*feval('fz',t(j),x(j),y(j));

kx2=h*feval('fx',t(j)+h/2,x(j)+kx1/2,y(j)+ky1/2,z(j)+kz1/2);
ky2=h*feval('fy',t(j)+h/2,x(j)+kx1/2,y(j)+ky1/2,z(j)+kz1/2);
kz2=h*feval('fz',t(j)+h/2,x(j)+kx1/2,z(j)+kz1/2);

kx3=h*feval('fx',t(j)+h/2,x(j)+kx2/2,y(j)+ky2/2,z(j)+kz2/2);
ky3=h*feval('fy',t(j)+h/2,x(j)+kx2/2,y(j)+ky2/2,z(j)+kz2/2);
kz3=h*feval('fz',t(j)+h/2,x(j)+kx2/2,z(j)+kz2/2);

kx4=h*feval('fx',t(j)+h,x(j)+kx3,y(j)+ky3,z(j)+kz3);
ky4=h*feval('fy',t(j)+h,x(j)+kx3,y(j)+ky3,z(j)+kz3);
kz4=h*feval('fz',t(j)+h,x(j)+kx3,z(j)+kz3);

x(j+1)=x(j)+(kx1+2*kx2+2*kx3+kx4)/6;
y(j+1)=y(j)+(ky1+2*ky2+2*ky3+ky4)/6;
z(j+1)=z(j)+(kz1+2*kz2+2*kz3+kz4)/6;
end

plot(t,x);

调用的fx,fy,fz分别为
function dx=fx(t,x,y,z)
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
dx=(1/2)*Gn*(z-nth)*x+(k/tin)*C1*cos(y-C2);
function dy=fy(t,x,y,z)
Gn=8.4e-13;n0=1.4e+24;tin=8e-12;k=0.015;a=6;
dy=(1/2)*a*Gn*(z-n0)-(k/tin)*(C1/x)*sin(y-C2);
function dz=fz(t,x,z)
Gn=8.4e-13;e=1.6e-19;V=1.0e-16;ts=2e-9;Ith=54e-3;n0=1.4e+24;
dz=(1.3*Ith)/(e*V)-(1/ts)*z-Gn*(z-n0)*x^2;

上边两个式子中有C1,C2,是两个数组，分别为

C1=[8.2e-011, -1.4679e-007, 2.2247e-007, 2.3087e-007, -1.0047e-007, 3.153e-008, 3.2047e-007, 1.0044e-007, -1.7736e-007, 1.8828e-007, -1.9584e-007, -1.1957e-007, 2.6189e-007, 2.4234e-008, 1.6046e-007, -1.2331e-007, -1.9984e-007, 4.6585e-007, -2.9163e-007, 1.4481e-008, -1.2686e-007, 6.6135e-007, -1.1262e-007, 1.3977e-007, -6.5892e-008, -3.6625e-007, -2.467e-007, 7.4485e-008, 1.475e-007, 2.7209e-007, 2.4435e-007, 1.0191e-007, 8.4376e-008, 1.5446e-008, -1.4789e-007, -1.9714e-007, 3.0642e-007, -3.1694e-007, 1.3063e-007, 1.2397e-007, -5.0877e-008, 8.2291e-008, -4.9396e-008, 4.4198e-008, -6.6902e-010, -2.0217e-008, 1.4079e-008, 2.1972e-008, 7.8652e-009, 1.0089e-008, 2.0888e-008, -5.0502e-008, 1.75e-007, 8.8072e-007, 1.6731e-006, 2.1356e-006, 2.0616e-006, 1.4339e-006, 2.6416e-007, -1.1631e-006, -2.3767e-006, -2.8328e-006, -1.3912e-006, 3.3939e-006, 1.7953e-005, 0.00010396, 0.00069155, 0.0052414, 0.045308, 0.44673, 5.0174, 64.156, 933.44, 15451, 2.9058e+005, 6.2072e+006, 1.5046e+008, 4.1381e+009, 1.2671e+011, 4.4426e+011, 8.9674e+010, 1.8253e+010, 4.2933e+009, 1.1719e+009, 3.7112e+008, 1.3615e+008, 5.7867e+007, 2.8467e+007, 1.6195e+007, 1.0647e+007, 8.0882e+006, 7.0847e+006, 7.1595e+006, 8.3524e+006, 1.1181e+007, 1.7265e+007, 3.0646e+007, 6.2556e+007, 1.4634e+008, 3.9367e+008, 1.2152e+009, 4.2877e+009, 1.7323e+010, 7.8254e+010, 2.8099e+011, 2.251e+011, 8.7931e+010, 3.4687e+010, 1.5524e+010, 8.0087e+009, 4.7717e+009, 3.2836e+009, 2.609e+009, 2.391e+009, 2.5263e+009, 3.0734e+009, 4.3131e+009, 6.934e+009, 1.2876e+010, 2.7086e+010, 6.4164e+010,
];
C2=[0, -5.3281, -9.5561, -12.689, -14.733, -15.694, -15.576, -14.385, -12.126, -8.8054, -4.4275, 1.0022, 7.4783, 14.996, 23.549, 33.133, 43.744, 55.374, 68.021, 81.677, 96.34, 112, 128.66, 146.31, 164.94, 184.56, 205.15, 226.71, 249.24, 272.73, 297.17, 322.57, 348.91, 376.2, 404.42, 433.58, 463.66, 494.67, 526.59, 559.44, 593.19, 627.84, 663.4, 699.85, 737.2, 775.43, 814.54, 854.54, 895.4, 937.14, 979.74, 1023.2, 1067.5, 1112.7, 1158.7, 1205.6, 1253.3, 1301.8, 1351.2, 1401.4, 1452.4, 1504.3, 1556.9, 1610.4, 1664.7, 1719.7, 1775.6, 1832.3, 1889.8, 1948, 2007, 2066.8, 2127.4, 2188.8, 2250.9, 2313.8, 2377.5, 2441.9, 2506.9, 2558.9, 2593.8, 2628.8, 2664.6, 2701.4, 2739, 2777.5, 2816.9, 2857.1, 2898.3, 2940.3, 2983.2, 3026.9, 3071.5, 3116.9, 3163.2, 3210.3, 3258.3, 3307.1, 3356.7, 3407.1, 3458.4, 3510.5, 3563.4, 3617, 3669.2, 3712.4, 3751.3, 3790.2, 3829.9, 3870.5, 3911.9, 3954.2, 3997.3, 4041.3, 4086.2, 4131.9, 4178.4, 4225.8, 4274, 4323, 4372.7,
];
我的意思就是想在主程序里，那个for 循环，每循环一次，就引用数组里的一个数值,主程序里算出来的数再接着继续存储在这两个数组里以便能继续算下去。拜托哪位高手帮帮我，我已经请教了很多人了，可是没弄出来，希望能有人帮忙，如果还有不清楚地，可以回个帖子，或者加我QQ也行，我的QQ：27866661，很着急，希望斑竹帮帮忙，推荐一下

你的t向量是0阿!调用的fx,fy,fz里似乎没有t阿
另外C1,C2是怎么算出来的
3000次循环而他们的维数只有121，是没算完吗？

如果每循环一次，就引用数组里的一个数值,主程序里算出来的数再接着继续存储在这两个数组里以便能继续算下去
那么可以把fx,fy,fz添加一个输入参量C

C1C2是另外的程序已经算出来的，怎么添加输入参量啊，这样说我也不太明白，能在我那个程序上直接改吗？谢谢啦！

调用的fx,fy,fz里似乎没有t阿
function dx=fx(t,x,y,z)
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
dx=(1/2)*Gn*(z-nth)*x+(k/tin)*C1*cos(y-C2);
c1,c2是数组算出来的dx也是数组，我不太明白你的具体问题
可以这样添加参量
function dx=fx(t,x,y,z，C1,C2)
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
dx=(1/2)*Gn*(z-nth)*x+(k/tin)*C1*cos(y-C2);
在主程序中
for j=1:M

kx1=h*feval('fx',t(j),x(j),y(j),z(j)，C1(j),C2(j));
.
.
.

我自己编的那个可能有点毛病，这样吧，你能帮我看看能不能编个程序，解这个微分方程组
dx(t)/dt=(1/2)*Gn*(z(t)-nth)*x(t)+(k/tin)*x(t-tt)*cos(y(t)-y(t-tt));
dy(t)/dt=(1/2)*a*Gn*(z(t)-n0)-(k/tin)*(x(t-tt)/x(t))*sin(y(t)-y(t-tt));
dz(t)/dt=(1.3*Ith)/(e*V)-(1/ts)*z(t)-Gn*(z(t)-n0)*x(t)^2;
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
n0=1.4e+24;k=0.015;a=6;
e=1.6e-19;V=1.0e-16;ts=2e-9;Ith=54e-3;tt=1.2e-9;
我只知道 用4阶龙格库塔法好像能解，别的我也不知道了，拜托斑竹一定帮我看看啊，已经找了好多人了，都没弄明白，郁闷阿

我是菜鸟，求教大家，如何实现将VB中的数值赋值给MATLAB中的函数？我的联系方式
QQ：409104373 E-mail:dongyunfeng2004@126.com

不知道算得对不对


function dy=fy(x,y,z,C1,C2)
Gn=8.4e-13;n0=1.4e+24;tin=8e-12;k=0.015;a=6;
dy=(1/2)*a*Gn*(z-n0)-(k/tin)*(C1/x)*sin(y-C2);


function dx=fx(x,y,z,C1,C2)
Gn=8.4e-13;nth=2.0e+24;k=0.015;tin=8e-12;
dx=(1/2)*Gn*(z-nth)*x+(k/tin)*C1*cos(y-C2);



function dz=fz(x,z)
Gn=8.4e-13;e=1.6e-19;V=1.0e-16;ts=2e-9;Ith=54e-3;n0=1.4e+24;
dz=(1.3*Ith)/(e*V)-(1/ts)*z-Gn*(z-n0)*x^2;



function [T,X,Y,Z]=qwe()




C1=[8.2e-011, -1.4679e-007, 2.2247e-007, 2.3087e-007, -1.0047e-007, 3.153e-008, 3.2047e-007, 1.0044e-007, -1.7736e-007, 1.8828e-007, -1.9584e-007, -1.1957e-007, 2.6189e-007, 2.4234e-008, 1.6046e-007, -1.2331e-007, -1.9984e-007, 4.6585e-007, -2.9163e-007, 1.4481e-008, -1.2686e-007, 6.6135e-007, -1.1262e-007, 1.3977e-007, -6.5892e-008, -3.6625e-007, -2.467e-007, 7.4485e-008, 1.475e-007, 2.7209e-007, 2.4435e-007, 1.0191e-007, 8.4376e-008, 1.5446e-008, -1.4789e-007, -1.9714e-007, 3.0642e-007, -3.1694e-007, 1.3063e-007, 1.2397e-007, -5.0877e-008, 8.2291e-008, -4.9396e-008, 4.4198e-008, -6.6902e-010, -2.0217e-008, 1.4079e-008, 2.1972e-008, 7.8652e-009, 1.0089e-008, 2.0888e-008, -5.0502e-008, 1.75e-007, 8.8072e-007, 1.6731e-006, 2.1356e-006, 2.0616e-006, 1.4339e-006, 2.6416e-007, -1.1631e-006, -2.3767e-006, -2.8328e-006, -1.3912e-006, 3.3939e-006, 1.7953e-005, 0.00010396, 0.00069155, 0.0052414, 0.045308, 0.44673, 5.0174, 64.156, 933.44, 15451, 2.9058e+005, 6.2072e+006, 1.5046e+008, 4.1381e+009, 1.2671e+011, 4.4426e+011, 8.9674e+010, 1.8253e+010, 4.2933e+009, 1.1719e+009, 3.7112e+008, 1.3615e+008, 5.7867e+007, 2.8467e+007, 1.6195e+007, 1.0647e+007, 8.0882e+006, 7.0847e+006, 7.1595e+006, 8.3524e+006, 1.1181e+007, 1.7265e+007, 3.0646e+007, 6.2556e+007, 1.4634e+008, 3.9367e+008, 1.2152e+009, 4.2877e+009, 1.7323e+010, 7.8254e+010, 2.8099e+011, 2.251e+011, 8.7931e+010, 3.4687e+010, 1.5524e+010, 8.0087e+009, 4.7717e+009, 3.2836e+009, 2.609e+009, 2.391e+009, 2.5263e+009, 3.0734e+009, 4.3131e+009, 6.934e+009, 1.2876e+010, 2.7086e+010, 6.4164e+010,
];
C2=[0, -5.3281, -9.5561, -12.689, -14.733, -15.694, -15.576, -14.385, -12.126, -8.8054, -4.4275, 1.0022, 7.4783, 14.996, 23.549, 33.133, 43.744, 55.374, 68.021, 81.677, 96.34, 112, 128.66, 146.31, 164.94, 184.56, 205.15, 226.71, 249.24, 272.73, 297.17, 322.57, 348.91, 376.2, 404.42, 433.58, 463.66, 494.67, 526.59, 559.44, 593.19, 627.84, 663.4, 699.85, 737.2, 775.43, 814.54, 854.54, 895.4, 937.14, 979.74, 1023.2, 1067.5, 1112.7, 1158.7, 1205.6, 1253.3, 1301.8, 1351.2, 1401.4, 1452.4, 1504.3, 1556.9, 1610.4, 1664.7, 1719.7, 1775.6, 1832.3, 1889.8, 1948, 2007, 2066.8, 2127.4, 2188.8, 2250.9, 2313.8, 2377.5, 2441.9, 2506.9, 2558.9, 2593.8, 2628.8, 2664.6, 2701.4, 2739, 2777.5, 2816.9, 2857.1, 2898.3, 2940.3, 2983.2, 3026.9, 3071.5, 3116.9, 3163.2, 3210.3, 3258.3, 3307.1, 3356.7, 3407.1, 3458.4, 3510.5, 3563.4, 3617, 3669.2, 3712.4, 3751.3, 3790.2, 3829.9, 3870.5, 3911.9, 3954.2, 3997.3, 4041.3, 4086.2, 4131.9, 4178.4, 4225.8, 4274, 4323, 4372.7,
];

a=0.000000012;
b=0.000000024;
M=3000;
h=(b-a)/M;
t=linspace(a,b,3001);
x=zeros(1,M+1);
y=zeros(1,M+1);
z=zeros(1,M+1);
t=a:h:b;
x(1)=6.4164e+010;
y(1)=4372.7;
z(1)=2.2129e+024;
for j=1:M

kx1=h*feval('fx',x(j),y(j),z(j),C1(j),C2(j));
ky1=h*feval('fy',x(j),y(j),z(j),C1(j),C2(j));
kz1=h*feval('fz',x(j),y(j));


kx2=h*feval('fx',x(j)+kx1/2,y(j)+ky1/2,z(j)+kz1/2,C1(j),C2(j));
ky2=h*feval('fy',x(j)+kx1/2,y(j)+ky1/2,z(j)+kz1/2,C1(j),C2(j));
kz2=h*feval('fz',x(j)+kx1/2,z(j)+kz1/2);


kx3=h*feval('fx',x(j)+kx2/2,y(j)+ky2/2,z(j)+kz2/2,C1(j),C2(j));
ky3=h*feval('fy',x(j)+kx2/2,y(j)+ky2/2,z(j)+kz2/2,C1(j),C2(j));
kz3=h*feval('fz',x(j)+kx2/2,z(j)+kz2/2);


kx4=h*feval('fx',x(j)+kx3,y(j)+ky3,z(j)+kz3,C1(j),C2(j));
ky4=h*feval('fy',x(j)+kx3,y(j)+ky3,z(j)+kz3,C1(j),C2(j));
kz4=h*feval('fz',x(j)+kx3,z(j)+kz3);


x(j+1)=x(j)+(kx1+2*kx2+2*kx3+kx4)/6;
y(j+1)=y(j)+(ky1+2*ky2+2*ky3+ky4)/6;
z(j+1)=z(j)+(kz1+2*kz2+2*kz3+kz4)/6;
C1=[C1 x(j+1)];
C2=[C2 y(j+1)];

end
T=t;
X=x;
Y=y;
Z=z;

plot(t,x);