Why is $1/r^2$ force law giving spiral trajectory?

Question

I have written a program to solve for Newton's 2nd Law of motion for a given force law, in 2D polar coordinates. It is known that if the force law is of the form $k/r^2$,we get conic sections as solution(https://en.wikipedia.org/wiki/Kepler_problem).

But when I run the program, I get a spiral trajectory (which isn't one of the conic sections).

Componentwise, Newton's 2nd law for central force $F(r)=k/r^2$(and unit mass) are-$$\frac{d^2r}{dt^2}=k/r^2$$ $$\frac{d^2\theta}{dt^2}=0$$ where the first equation is the radial part and the 2nd equation is for angular part(which vanishes for Central force problems,due to conservation of angular momentum).

The equations to be used in the program are a set of 4 coupled first order ODEs, derived from the above 2nd order ODEs.(Since odeint cant directly handle 2nd order ODE).The equations I used in the program were--- $$dr/dt=v_r$$$$dv_r/dt=k/r^2$$$$d\theta/dt=v_\theta$$$$dv_\theta/dt=0$$with $k=1$ for simplicity.

I have tried with many possible initial conditions, and tried running for different values of $k$ to ensure. Why is this happening? Is this because of the choice of initial condition or some faulty implementation or some numerical problem? Or is their a problem with the physics?

Below is my minimal code

def vec(w, t):
    r, vr, theta, vtheta = w
    return [vr, r**(-2.0), vtheta, 0]

def newton(vec, initial, t):
    wsol = odeint(vec, initial, t)
    return [wsol[:, 0], wsol[:, 2]]

T = np.linspace(0, 50, 1000)    
initial = [2, 10, radians(0),2] #The order is radius, radial velocity, angle, angular velocity
R = newton(vec, initial, T)[0]
Theta = newton(vec, initial, T)[1]
plt.polar(Theta, R, "r", lw="1")
plt.show()

Answer 1

The correct dynamic equations in the polar coordinates should be

$ \dot{v_r} = \omega^2 r - \alpha/r^2 \\ \dot{\omega} = - 2 v_r \omega /r\\ \dot{\theta} = \omega \\ \dot{r} = v_r $

Here is the fixed Python code:

from math import *
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint


def vec(w,t):
    r,vr,theta,omega=w
    return [vr, r*omega**2 - 1e0/r**2, omega, -2*vr*omega/r]

def newton(vec,initial,t):
    wsol=odeint(vec,initial,t)
    return [wsol[:,0],wsol[:,1],wsol[:,2],wsol[:,3]]

t=np.linspace(0,4,1121)
#The order is radius,radial velocity,angle,angular velocity.
initial=[0.5, -1e0, radians(0), 1e0] 

r=newton(vec,initial,t)[0]
vr=newton(vec,initial,t)[1]
theta=newton(vec,initial,t)[2]
omega=newton(vec,initial,t)[3]

plt.polar(theta,r, "r", lw="1")
plt.show()

Here is the output: