function a = lrp(t,r)

global po sigma1 Ro1 k amp mul f rhol Gs mus

%w=sin(2.0*pi*(t)*f)+sin(2.0*7*pi*(t)*f);
w=sin(2.0*pi*(t)*f);

a = [r(2);...
    (...
    (...
    (po+2.0*sigma1/Ro1)*Ro1^(3*k)*r(1)^(-3*k) ...
    -po ...
    -amp*w ...
    -Gs*(r(1)-Ro1)*r(1)^(-4)...
    -mus*r(2)*r(1)^(-4)...
    -2.0*sigma1*r(1)^(-1)...
    -4.0*r(2)*mul*r(1)^(-1)...
 )/rhol ...
    -1.5*r(2)^2 ...
    )*r(1)^(-1) ...
];