请看一下方程的问题

zhengdong001

function dx=newvdp(t,x)
T1=20;            
I1=1.3696e-2;            
I2=2.92645e-4;
I3=6.7308e-3;
I4=3.5117e-3;
I5=3.0142e-1;
I6=1.0272;
c1=0.6233;            
c2=18.412;
c3=21.65;
c4=22.55;
c5=15.394;
k1=1.356e4;            
k2=9.975e8;
k3=1.017e6;
k4=1.469e6;
k5=1.4875e9;
R1=40;            
R2=160;
R3=35;
R4=180;
em1=1;         
em2=1;
e1=em1+3.45*cos(20*t);   
e2=em2+0.26*cos(30*t);
w1=c2*(R1*x(4)-R2*x(6)-diff('e1','t',1))+k2*(R1*x(3)-R2*x(5)-e1);   
w2=c4*(R3*x(8)-R4*x(7)-diff('e2','t',1))+k4*(R3*x(7)-R4*x(9)-e2);
dx=zeros(12,1);   
dx(1)=x(2);         
dx(2)=(T1-c1*(x(2)-x(4))-k1*(x(1)-x(2)))/I1;
dx(3)=x(4);
dx(4)=(-c1*(x(4)-x(2))-k1*(x(2)-x(1))-R1*w1)/I2;
dx(5)=x(6);
dx(6)=(R2*w1-c3*(x(6)-x(8))-k3*(x(5)-x(7)))/I3;
dx(7)=x(8);
dx(8)=(-c3*(x(8)-x(6))-k3*(x(7)-x(5))-R3*w2)/I4;
dx(9)=x(10);
dx(10)=(R4*w2-c5*(x(10)-x(12))-k5*(x(9)-x(11)))/I5;
dx(11)=x(12);
dx(12)=(-c5*(x(12)-x(10))-k5*(x(11)-x(9)))/I6;


[T,x]=ode45('newvdp',[0,0.5],[0,0,0,0,0,0,0,0,0,0,0,0]);
plot(T,x(:,1),'r',T,x(:,2),'b*',T,x(:,3),'k-.')

大家看一下这个方法有什么问题，我算了十多个小时也没结束运算，目前还在算？