帮个忙吧，用ode45解微分方程组，刚开始学，请指点一二！！

西西想幸福

这是我的方程程序：
function dx=modfun(t,x);
global rc d1 r1 d2 r2 cpg gms cps h1 h2
rc=rc(x(3));
d1=d1(x(4));
r1=r1(x(1),x(3));
d2=(x(4));
r2=r2(x(2),x(3));
cpg=cpg(x(4));
gms=gms(x(3));
cps=cps(x(5));
h1=h1(x(5));
h2=h2(x(5));
dx(1)=-0.99/1.98*10^8*r1 ;                                                             %定义第一个方程，根据反应中氢气的质量守恒得出
dx(2)=-0.99/1.98*10^8*r2;                                                               %定义第二个方程，根据反应中一氧化碳的质量守恒得出
dx(3)=-0.99/4.56*10^4*[r1+r2];                                                          %根据固相反应量的出的第三个方程，根据固相的质量守恒得出
dx(4)=14.07*(x(5)-x(4))/cpg);                                                      %为计算方便，将气相温度tg换为x4，后面也是如此（化简后）,根据能量守恒得出方程，并未考虑热量损失
dx(5)=0.99/gms/cps*[pi*4*10^4*-h1*r1-h2*r2];                                     %固相温度tsol也换为x5，根据固相的能量守恒
function rc(x3)
rc=(0.125-3.78*x3)^1/3
function d1(x4)
d1=1.467*10^-6*x4^1.75
function r1(x1,x3)
r1=-[4*pi*rc^2*(0.529-x1)]/[1+k1+rc(x3)/d1-4*rc(x3)/d1]

function d2(x4)
d2=1.467*10^-6*x4^1.75  

function r2(x2,x3)
r2=-[4*pi*rc^2*(0.347-x2)]/[1+k2+rc(x3)/d2-4*rc(x3)/d2]
   
function cpg(x4)
cpg=7.045+0.0017*x4-2.603*10^4*x4^-2
function gms(x3)
gms=2527.55*[692-447*rc(x3)^3]/[692-341*rc(x3)^3]
function cps(x5)
cps=[0.6975+1.8058*x5-8.462*x5*rc(x3)^3+102.234*rc(x3)^3]/[0.3875-0.3*rc(x3)^3]
function h1(x5)
h1=80342.124-8.56*x5+2.1*10^-3*x5^2+1.2*10^4*x5^-1
function h2(x5)
h2=-8621.38-5.56*x5+1.043*10^-3*x5^2-1.98*10^5*x5^-1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
tspan=[0,1000];                                                                          %定义变量求解区间
x0=[1.4e-6,1.5e-6,0,650,298];                                                            %设置初值，并设在z=0处x1的值是1.4*10^-6，x2的值是1.5*10^-6,tg是650k
options=odeset('RelTol',1e-6);                                                           %设置相对误差
[t,x]=ode45(@modfun,tspan,x0,options);                                                  %调用ode45函数解微分方程
figure;
plot(t,x(:,1),'k-');
hold on;
plot(t,rc,'k:');
set(gca,'Fontsize',12);
xlabel('\itt','Fontsize',16);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
我只是照着书上和一些帖子里面的程序编的，有好多东西都不太懂！！说出现了错误我也不知道怎么改！！请各位指点一下吧！！谢谢了！！哦，对了！！我的错误代码是：
??? Error using ==> feval
Undefined function or method 'modfun' for input arguments of type 'double'.
Error in ==> odearguments at 111
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in ==> modfun1 at 6
[t,x]=ode45(@modfun,tspan,x0,options);                                                  %调用ode45函数解微分方程



帮帮忙吧！！十分感谢！！

西西想幸福

都没人理我啊！！~~！！

321forever

Undefined function or method 'modfun' for input arguments of type 'double'.
每个function 都要保存，且存在同一个文件夹中

西西想幸福

回复 3 # 321forever 的帖子

我已经保存了啊！！可是还是不行啊？是不是这个方法不能解这个方程组啊？那要用什么方法呢？

西西想幸福

回复 4 # 西西想幸福 的帖子

帮帮忙看一下吧，各位！！先谢过了！！

321forever

function dx=modfun(t,x);
去掉；
dx(4)=14.07*(x(5)-x(4))/cpg);   多个括号

西西想幸福

回复 6 # 321forever 的帖子

改过了，可是错误还是这个!!谢谢了！！再帮我看下？真的麻烦你了！！


??? Error using ==> feval
Undefined function or method 'modfun' for input arguments of type 'double'.

Error in ==> odearguments at 111
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> modfun1 at 6
[t,x]=ode45(@modfun,tspan,x0,options);                                                  %调用ode45函数解微分方程

321forever

1. lz把所有的[]改成（）
2. k1和k2 是什么？
3. 在modfun 中加入dx=zeros(5,1)，或是在最后加上dx=dx(:) 还有return
4. 在plot(t,rc,'k:'); rc在x中是第几个，假如是第3个，plot(t,x(:,3),'k:');
都改过了应该可以出图
Error using ==> feval
Undefined function or method 'modfun' for input arguments of type 'double'.

这个问题我试了下没有出现 lz再试试吧
function dx=modfun(t,x)
rc=(0.125-3.78*x(3))^1/3;
d1=1.467*10^-6*x(4)^1.75;
r1=-(4*pi*rc^2*(0.529-x(1)))/(1+k1+rc/d1-4*rc)/d1;
d2=1.467*10^-6*x(4)^1.75;
r2=-(4*pi*rc^2*(0.347-x(2)))/(1+k2+rc/d2-4*rc/d2);
cpg=7.045+0.0017*x(4)-2.603*10^4*x(4)^-2;
gms=2527.55*(692-447*rc^3)/(692-341*rc^3);
cps=(0.6975+1.8058*x(5)-8.462*x(5)*rc^3+102.234*rc^3)/(0.3875-0.3*rc^3);
h1=80342.124-8.56*x(5)+2.1*10^-3*x(5)^2+1.2*10^4*x(5)^-1;
h2=-8621.38-5.56*x(5)+1.043*10^-3*x(5)^2-1.98*10^5*x(5)^-1;
dx=zeros(5,1);
dx(1)=-0.99/1.98*10^8*r1 ;                   %定义第一个方程，根据反应中氢气的质量守恒得出
dx(2)=-0.99/1.98*10^8*r2;                    %定义第二个方程，根据反应中一氧化碳的质量守恒得出
dx(3)=-0.99/4.56*10^4*(r1+r2);               %根据固相反应量的出的第三个方程，根据固相的质量守恒得出
dx(4)=14.07*(x(5)-x(4))/cpg;                 %为计算方便，将气相温度tg换为x4，后面也是如此（化简后）,根据能量守恒得出方程，并未考虑热量损失
dx(5)=0.99/gms/cps*(pi*4*10^4*-h1*r1-h2*r2); %固相温度tsol也换为x5，根据固相的能量守恒
return

这是改的lz的程序，删了其他的function，只留了这个modfun，k1和k2是什么？

西西想幸福

不好意思，昨晚有事就先走了！！给我指出这么多错误，我试下！！嘻嘻！！k1和k2 也是两个参数，就跟d1和d2 差不多！！以前改的时候不小心弄没了！！我加上试下！！真的很感激！！

西西想幸福

好像可以啊！！就是好慢，还在算！不过没有红色字体的报错信息，我已经开心的要死了！！真的真的很感谢！！

西西想幸福

为什么已经算了快要一个小时了，怎么还没有算完？是不是我的方程有问题？还是因为t是[0,1000]所以算的时间比较长？有没有人知道这个？在帮帮忙吧！！

321forever

回复 11 # 西西想幸福 的帖子

lz能把改好的加上k1和k2的程序再发下么，还有原方程发上来看看

西西想幸福

我又把所有的参数从新算了一遍，更正了一下！！然后再算！！刚算还不知道结果怎么样@！！恩，程序如下：
clc
clear
tspan=[0,1000];                                                                          %定义变量求解区间
x0=[1.4*10^-6,1.5*10^-6,0,650,298];                                                            %设置初值，并设在z=0处x1的值是1.4*10^-6，x2的值是1.5*10^-6,tg是650k
options=odeset('RelTol',1e-6);                                                           %设置相对误差
[t,x]=ode45(@modfun,tspan,x0,options);                                                  %调用ode45函数解微分方程
figure;
plot(t,x(:,1),'k:');
hold on;
plot(t,x(:,3),'k:');
set(gca,'Fontsize',12);
xlabel('\itt','Fontsize',16);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dx=modfun(t,x)
rc=(0.125-3.78*x(3))^1/3;
d1=1.467*10^-6*x(4)^1.75;
d2=3.828*10^-7*x(4)^1.75;
k1=0.225*exp(-179.14/x(4));
k2=0.650*exp(-342.4324/x(4));
r1=-(4*pi*rc^2*(0.529-x(1)))/(1/k1+rc/d1-4*rc^2/d1);
r2=-(4*pi*rc^2*(0.347-x(2)))/(1/k2+rc/d2-4*rc^2/d2);
cpg=104.1273+0.9806*x(4)-363688*x(4)^-2-314.5*x(1)-655.2*x(2)-0.3224*x(1)*x(4)-2.8364*x(4)*x(2)+24000*x(4)^-2*x(1)+9372000*x(4)^-2*x(2);
gms=2533.3*(-4.4864*rc^3+0.0069)/(-0.0342*rc^3+0.0069);
cps=(3.7125+0.006*x(5)*rc^3+79.9140*rc^3)/(0.3875-0.3*rc^3);
h1=999.8621-3.79*x(5)+0.0011*x(5)^2+1.2*10^4*x(5)^-1;
h2=840.3932-0.79*x(5)+5.333*10^-4*x(5)^2-1.95*10^5*x(5)^-1
dx(1)=-0.99/1.98*10^8*r1 ;                                                             %定义第一个方程，根据反应中氢气的质量守恒得出
dx(2)=-0.99/1.98*10^8*r2;                                                               %定义第二个方程，根据反应中一氧化碳的质量守恒得出
dx(3)=-0.99/4.56*10^4*(r1+r2);                                                          %根据固相反应量的出的第三个方程，根据固相的质量守恒得出
dx(4)=14.07*(x(5)-x(4))/cpg;                                                           %为计算方便，将气相温度tg换为x4，后面也是如此（化简后）,根据能量守恒得出方程，并未考虑热量损失
dx(5)=0.99/gms/cps*(pi*4*10^4*(x(5)-x(4))-h1*r1-h2*r2);                                %固相温度tsol也换为x5，根据固相的能量守恒
dx=dx(:);
return
不知道是不是t的取值范围有点大，现在我的conmand窗口里面一直都在显示h2=。。。。而且一直都在变化，不知道这是不是正确的现象！！嘻嘻，真的很感激你啊！！

321forever

回复 13 # 西西想幸福 的帖子

客气啦，h2方程后面少了个；刚刚试了下，也没有运行出来。问下哈，lz能贴下原方程么。

西西想幸福

原方程就是一组动力学微分方程：
u*dx1dz+np*r1=0      
u*dx2dz+np*r2=0   
dTgdz-(np*Ap*h*(Ts-Tg))/Gmg/Cpg=0
u*dx3dz+np*(r1+r2)=0
dTsdz-np/Gms/Cps*（Ap*h*(Ts-Tg)-h1*r1-h2*r2）=0
其中u、np、Ap、Gmg和r1、r2是常数
不知道是不是要这个原方程！！我把这些常数带入后，得出我编的那个方程！然后还有那些关系式！！这样的方程ode45能解么？