急急急！请教前辈，高手，怎么求极限环啊

冰雪之巅

下面是我自己编写的M文件，我用ODE45求它的相轨迹怎么会逐渐增大呢？是哪里的问题啊？能用ODE45求极限环吗？如果不能极限环怎么来求？？？请前辈，高手指教！！！！
function  xdot=mocali(t,x)
xdot=zeros(4,1);
syms k1 k2 k3 k4 k5 c1 c2 c3 m1 m2;
k1=300000;k2=200000;k3=500000;k4=1000000;k5=1000000;a=(pi/2-0.2);
b=(pi/2-0.2);u=0.4;n=50;
c1=1;c2=3;c3=3;m1=1;m2=10;
xdot(1)=x(2);
xdot(2)=-(k1+k2)*x(1)/m2-(c1+c2)*x(2)/m2+k1*x(3)/m2+c1*x(4)/m2-k4*(x(1)-x(3))^2/m2-k5*(x(1)-x(3))^3/m2+u*n*sin(pi/2-0.2)/m2;
xdot(3)=x(4);
xdot(4)=-k1*x(1)/m1-c1*x(2)/m1-(k3-k1)*x(3)/m1-(c3-c1)*x(4)/m1-k4*(x(1)-x(3))^2/m1-k5*(x(1)-x(3))^3/m1+(n*cos(pi/2-0.2)-u*n*sin(pi/2-0.2))/m1;


>>t0=0;tf=20;x0=[-0.0008 0.005 -0.0004 0.001]';
>>[t,x]=ode45('mocali',[t0,tf],x0);