matlab数值求解边界条件的微分方程组

matlab数值求解边界条件的微分方程组
解一个方程组,文献上给出了如下的形式：
dx/dz=-a1*x-g1*v1/v2*x(m+n+2d)+c1*y
dy/dz=a1*y-g1*v1/v2*x(m+n+2d)-c1*x
dm/dz=-a2*m+g1*(m+d)*(x+y)-g2*v2/vs*m*(w+v)+c2*n
dn/dz=a2*n-g1*(n+d)*(x+y)+g2*v2/vs*n*(w+v)-c2*m
dw/dz=as*v-g2*v2/vs*v*(m+n)-gb*w*v
dv/dz=-as*w+g2*v2/vs*w*(m+n)-gb*w*v
其中z是自变量,x y m n w v 是函数,其他的都是常数.
边界条件和初值条件为：
x(0)=0.47 y(50)=0.47 w(0)=0.0075 v(0)=0.7e-6
不知道4个条件能不能做.
按照matlab求解边界条件微分方程bvp4c方法,我进行了如下处理：
y(1)=x,y(2)=y,y(3)=m,y(4)=n,y(5)=w,y(6)=v
方程M文件
function dydx=ivode(x,y)
dydx=[-a1*y(1)-g1*v1/v2*y(1)*(y(3)+y(4)+2*d)+c1*y(2);
a1*y(2)+g1*v1/v2*y(2)*(y(3)+y(4)+2*d)-c1*y(1);
-a2*y(3)+g1*(y(3)+d)*(y(1)+y(2))-g2*v2/vs*y(3)*(y(5)+y(6))+c2*y(4);
a2*y(4)-g1*(y(4)+d)*(y(1)+y(2))+g2*v2/vs*y(4)*(y(5)+y(6))-c2*y(3);
as*y(6)-g2*v2/vs*y(6)*(y(3)+y(4))-gb*y(5)*y(6);
-as*y(5)+g2*v2/vs*y(5)*(y(3)+y(4))-gb*y(5)*y(6)];
边界条件x(0)=0.47 y(50)=0.47 w(0)=0.0075 v(50)=0.7e-6
边界条件M文件
function res=ivbc(ya,yb)
res=[ya(1)-0.47;ya(5)-0.0075;yb(2)-0.47;yb(6)-0.7e-6];
各种常数参数如下：
a1=0.31;a2=0.25;as=0.2;
g1=0.51;g2=0.38;gb=5e-14;
h=6.62e-34;kb=1.38e-23;t=298;
c1=1.0e-5;c2=6.0e-5;
v1=3e8/1365e-9;v2=3e8/1455e-9;vs=3e8/1550e-9;
dv2=0.4e-12;
d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)*kb*t)-1));
如果需要6个边界条件的话,暂定ya(3)=0,yb(4)=0.最好有完整的程序,自己调了半天全是错.
1楼的大神指出的错误是我自己没打对,
d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)/(kb*t))-1));
算出d的值为1.2298e-31.

请核实一下条件：
d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)*kb*t)-1));
这个式子确定无误吗?
按照现在给的数值：h=6.62e-34;kb=1.38e-23;
h*(v1-v2)*kb*t的值约为3.7e-041,所以exp(h*(v1-v2)*kb*t)-1=0,d为inf,根本没法往下算.
再问： 感谢大神指出错误，是我自己打错了，现在更正过来， d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)/(kb*t))-1)); 希望大神指点一下接下怎么做，看了不少MATLAB的书，很少有教怎么解边界条件下微分方程的，求指点啊。
再答： function zd513599949 % 解的初始估计 solinit = bvpinit(linspace(0,50,5),[1 0 0 0 0 0]); % BVP问题求解 sol = bvp4c(@ivode,@ivbc,solinit); % 结果绘图 t=sol.x; vars={'x', 'y', 'm', 'n', 'w', 'v'}; for i=1:length(vars) subplot(2,3,i); plot(t,sol.y(i,:)); xlabel('t'); ylabel(vars{i}); end function dydx=ivode(x,y) % 常数定义 a1=0.31;a2=0.25;as=0.2; g1=0.51;g2=0.38;gb=5e-14; h=6.62e-34;kb=1.38e-23;t=298; c1=1.0e-5;c2=6.0e-5; v1=3e8/1365e-9;v2=3e8/1455e-9;vs=3e8/1550e-9; dv2=0.4e-12; % d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)*kb*t)-1)) d=2*h*v2*dv2*(1+1/(exp(h*(v1-v2)/(kb*t))-1)); % 微分方程 dydx=[-a1*y(1)-g1*v1/v2*y(1)*(y(3)+y(4)+2*d)+c1*y(2); a1*y(2)+g1*v1/v2*y(2)*(y(3)+y(4)+2*d)-c1*y(1); -a2*y(3)+g1*(y(3)+d)*(y(1)+y(2))-g2*v2/vs*y(3)*(y(5)+y(6))+c2*y(4); a2*y(4)-g1*(y(4)+d)*(y(1)+y(2))+g2*v2/vs*y(4)*(y(5)+y(6))-c2*y(3); as*y(6)-g2*v2/vs*y(6)*(y(3)+y(4))-gb*y(5)*y(6); -as*y(5)+g2*v2/vs*y(5)*(y(3)+y(4))-gb*y(5)*y(6)]; % 边界条件 function res=ivbc(ya,yb) res=[ya(1)-0.47; ya(5)-0.0075; yb(2)-0.47; ya(6)-0.7e-6; ya(3); yb(4)];