微分方程里嵌套quadl到底行不行?

sundayfen

各位大侠，帮忙看看这程序解不解的出来，我快疯了，好几次都是电脑运行了好长时间，没出结果就瘫痪了
clear all; clc
[t,y]=ode45(@stokes,[0,1e1],[1e-4,0,1e-1]);
plot(t,y(:,1),'r',t,y(:,2),'b')


function dy=stokes(t,y)
h=6.63*1e-34; gamma=0.71; c=3e8; ls=15e-3; l=ls; n1=2.183;
R1L=0.999; R2L=0.998; R1s=0.928; R2s=0.97;
lc=(n1*l+(0.75-l)); tr=2*lc/c;
tL=tr/(0.05-log(R1L*R2L)); ts=tr/(0.05-log(R1s*R2s));
sigma=2.5e-22; ksp=2e-6;

wp=300e-6; wr=300e-6; wl=300e-6; ws=300e-6;
krs=wr^2/ws^2; klr=wl^2/wr^2; tao1=98e-6; g=4.5*1e-11;
mius=c*(10-1.064*0.89)/(1.064e-5); G=2*g*h*mius*c*ls; n0=1e-4;

dy=[4.*sigma.*l./wl.^2./tr.*y(1).*n0.*quadl(@(r)(exp(-2.*r.^2.*(1/wp^2+1/wl^2)).*exp(-gamma.*sigma.*c.*exp(-2.*r.^2./wl.^2).*y(3)).*2.*r),0,1e3)-1./(tr.*(1+krs)).*G.*y(1).*y(2)-y(1)./tL;krs./(tr.*(1+krs)).*klr.*G.*y(1).*y(2)-y(2)./ts;y(1)];

我比较怀疑的是这句 quadl(@(r)(exp(-2.*r.^2.*(1/wp^2+1/wl^2)).*exp(-gamma.*sigma.*c.*exp(-2.*r.^2./wl.^2).*y(3)).*2.*r),0,1e3)，积分区间是本来是0到无穷， 但是我怕quadl没法用，就差不多估计了个数。但是，算不出来……  哪位高人帮忙看看？

高人哪，帮忙看看吧
微分方程里嵌套quadl到底行不行啊

[ 本帖最后由 ChaChing 于 2010-6-13 18:42 编辑 ]