I am creating a dynamics question where an asteroid in space does not have enough tangential velocity to escape the pull of a planet, so it gradually will spiral around the planet eventually crashing into it. Now the question would ask the student to just find the position of the asteroid at a specific time $t = t_f\ \rm{s}$ where $t_f$ would be randomized such as the asteroid would not have crashed into the planet yet.
I started my derivation from Newton's Gravitational Law and the kinematics equation $v_r.dv_r=a_r.dr$, such as:
$${v_r.dv_r=\frac{-G.M}{r^2}.dr}$$
$$\int_0^{v_r}{v_r}.dv_r=-\int_{r_0}^r{\frac{G.M}{r^2}}$$
$$v_r^2=2.G.M.(\frac{1}{r}-\frac{1}{r_0})$$
Now, this is a non-linear first-order ode that I have tried everything I know to solve but failed miserably. MATLAB ODE solvers do not like this formulation as it does not contain the independent variable in the formulation. What ends up happening is it plugs the initial condition $r(0) = r_0 \rm \ m$, which gives back a velocity of $0\ \rm m/s$, which in turn would keep $r$ at $r_0$ repeating the cycle. I was wondering if there is some sort of solution, whether numerical or analytical to solve this kind of problem.
I am going to include the very basic code I wrote when trying to solve this problem.
R = 3*10^6;
d0 = 42000;
G= 6.6743*10^-11;
tspan= [0,1000];
r0= R+d0;
[t,r] = ode113(@(t,r) -(2*G*M*((1/r)-(1/(R+d0))))^0.5, tspan, r0);
plot(t,r);