function dr=rp(t,r)
global R0 sigema rou0 Pv gamma yita P0 f1 f2 P1 P2
dr=zeros(2,1);
dr(1)=r(2);
dr(2)=(1/r(1))*((((P0+2*sigema/R0-Pv)*(R0/r(1))^(3*gamma)+Pv-2*sigema/r(1)-4*yita*r(2)/r(1)-P0-P1*cos(2*pi*f1*t)-P2*cos(2*pi*f2*t)))/rou0-1.5*r(2)^2);