ode45求解微分方程组，一直报错，请各位大侠帮我！

姚明

（1）求解问题为：解8个方程构成的微分方程组，微分方程的系数是时变的。
     （2）我的解题思路是：在M文件建立函数“function dx=myfun8(t,x)”，再利用“[t,x]=ode45('myfun8',[0,0.02],[0 0 0 0 0 0 0 0])”求解，初值均为0。
     （3）程序问题：（a）用ode45显示结果为“out of memory”,并且用时明显很长；
                     (b)若荣ode15s等解刚性方程组的命令，显示结果为“Warning: Failure at t=7.619368e-006.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.706943e-020) at time t.
> In ode15s at 738”
                    （c）若改变初值，比如把后两个初值改为：0.1、0.1，可以用ode15s求出结果，但初值改变了，也不知道得到的结果是否准确。
     纠结了很久了，实在是初学没有经验，麻烦大侠多指导我，不胜感激。


程序代码为：
function dx=myfun8(t,x)
U1=[7.9677e3*sin(100*pi*t); 7.9677e3*sin(100*pi*t+1.5*pi); 7.9677e3*sin(100*pi*t-1.5*pi); 0;  0;  0;];
%定义初值
Lp1=1e-3; Lp2=1e-3; Lp3=1e-3;
rp1=2.5;    rp2=2.5;    rp3=2.5;   
rc1=2274;   rc2=2274;   rc3=2274;
Np1=78;    Np2=78;    Np3=78;
Lc1=3.017;      Lc2=3.017;     Lc3=3.017;     Lc4=3.160;       Lc5=3.160;       Lc6=3.160;       Lc7=3.160;
Ac1=0.5555;     Ac2=0.5555;    Ac3=0.5555;    Ac4=0.5635;      Ac5=0.5635;       Ac6=0.5635;     Ac7=0.5635;

Nv=[Np1 0 0; 0 Np2 0; 0 0 Np3;];
Nf=[Np1 -Np2 0; 0 -Np2 Np3;];
C=[1 0; -1 -1; 0 1;];

ud1=1/(eps+0.61425*cosh(4.55*x(7)/Ac1));  ud2=1/(eps+0.61425*cosh(4.55*x(8)/Ac2));  ud3=1/(eps+0.61425*cosh(4.55*(-x(7)-x(8))/Ac3));  
ud4=1/(eps+0.61425*cosh(4.55*x(7)/Ac4));  ud5=1/(eps+0.61425*cosh(4.55*x(7)/Ac5));  ud6=1/(eps+0.61425*cosh(4.55*(-x(7)-x(8))/Ac6));  ud7=1/(eps+0.61425*cosh(4.55*(-x(7)-x(8))/Ac7));

Rmd1=Lc1/(ud1*Ac1+eps);  Rmd2=Lc2/(ud2*Ac2+eps);  Rmd3=Lc3/(ud3*Ac3+eps);  Rmd4=Lc4/(ud4*Ac4+eps);  Rmd5=Lc5/(ud5*Ac5+eps);  Rmd6=Lc6/(ud6*Ac6+eps);  Rmd7=Lc7/(ud7*Ac7+eps);
a1=Rmd1*Rmd2+Rmd1*Rmd3+Rmd1*Rmd6+Rmd1*Rmd7+Rmd2*Rmd3+Rmd2*Rmd6+Rmd2*Rmd7+Rmd4*Rmd2+Rmd4*Rmd3+Rmd4*Rmd6+Rmd4*Rmd7+Rmd5*Rmd2+Rmd5*Rmd3+Rmd5*Rmd6+Rmd5*Rmd7;
Rn=[(Rmd2+Rmd3+Rmd6+Rmd7)/(a1+eps)  -Rmd2/(a1+eps); -Rmd2/(a1+eps)  (Rmd1+Rmd2+Rmd4+Rmd5)/(a1+eps);];
Ld=Nv*C*Rn*Nf+eps;

r=[rp1 0 0 0 0 0;  0 rp2 0 0 0 0;  0 0 rp3 0 0 0;  rc1 0 0 -rc1 0 0;  0 rc2 0 0 -rc2 0;  0 0 rc3 0 0 -rc3;];
a2=Ld(1,1)*Ld(2,2)*Ld(3,3)-Ld(1,1)*Ld(2,3)*Ld(3,2)-Ld(2,1)*Ld(1,2)*Ld(3,3)+Ld(2,1)*Ld(1,3)*Ld(3,2)+Ld(3,1)*Ld(1,2)*Ld(2,3)-Ld(3,1)*Ld(1,3)*Ld(2,2)+eps;
L=[1/Lp1,0,0,1/Lp1,0,0;  0,1/Lp2,0,0,1/Lp2,0;  0,0,1/Lp3,0,0,1/Lp3;  0,0,0,-(Ld(2,2)*Ld(3,3)-Ld(2,3)*Ld(3,2))/a2,(Ld(1,2)*Ld(3,3)-Ld(1,3)*Ld(3,2))/a2,(-Ld(1,2)*Ld(2,3)+Ld(1,3)*Ld(2,2))/a2;...
0,0,0,(Ld(2,1)*Ld(3,3)-Ld(2,3)*Ld(3,1))/a2,-(Ld(1,1)*Ld(3,3)-Ld(1,3)*Ld(3,1))/a2,(Ld(1,1)*Ld(2,3)-Ld(1,3)*Ld(2,1))/a2;  0,0,0,(-Ld(2,1)*Ld(3,2)+Ld(2,2)*Ld(3,1))/a2,(Ld(1,1)*Ld(3,2)-Ld(1,2)*Ld(3,1))/a2,-(Ld(1,1)*Ld(2,2)-Ld(1,2)*Ld(2,1))/a2;];

dx(1:6)=L*(U1-r*[x(1); x(2);  x(3); x(4);  x(5); x(6);]);
dx(7:8)=Rn*Nf*[dx(4); dx(5); dx(6);];
dx=dx(:);
end