Using scipy.integrate.ode to solve protoplanetary hydrostatic equilibria

Question

I have written code to solve the protoplanetary hydrostatic equilibria equations using first order euler method, and the results are as expected. I now want to rewrite it using fourth order Runge-Kutta solution, which is builtin to scipy.integrate.ode

Expected result, attained by using a custom written euler method:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode


universalGasConstant = 8.31
gravitationalConstant = 6.67 * 10**-11
astronomicalUnit = 1.5 * (10**11)
coreMass = 5.9 * 10**24
sunMass = 2 * (10**30)
coreRadius = 0
coreDensity = 5500.0

semiMajorAxis = 5 * astronomicalUnit

gasMeanMolecularMass = 2.3 * 10**-3

coreRadius = (coreMass / ((4 / 3) * np.pi * coreDensity))**(1 / 3)
print('coreradius is', coreRadius)

hillSphere = (1**-4)*(semiMajorAxis) * (((coreMass) / (3 * sunMass))**(1 / 3))
print('hillsphere is', hillSphere)#THIS IS THE HILLSPHERE

numberOfSlices = 100000000 #NUMBER OF SLICES

deltaRadius = hillSphere / numberOfSlices
print('slice length is', deltaRadius)

#gravitationalForce = ((gravitationalConstant * coreMass) / 
    #(coreRadius + (deltaRadius)))


sunSurfaceTemperature = 5800
sunRadius = 6.96 * 10**8
equilibriumTemperature = sunSurfaceTemperature * \
    (sunRadius / (2 * semiMajorAxis))**(1 / 2)
print('equilibrium temperature at', semiMajorAxis,
      'metres is', equilibriumTemperature, 'kelvin')



pressureAtRadius = 0
initialDensity = 10.0 #INITIAL DENSITY
densityAtRadius = 0

muOverRT = gasMeanMolecularMass/(equilibriumTemperature * universalGasConstant)
print('muOverRT is', muOverRT)

totalMass = coreMass

densityAtRadius = initialDensity

def function(y, r):
        totalMass = y[0]
        distanceFromCOG = y[1]

        # the differential equations
        dRhodr = -(((gravitationalConstant*totalMass*
                       densityAtRadius*muOverRT))/(distanceFromCOG**2))
        dMdr = 4*np.pi*gravitationalConstant*densityAtRadius
        #dPdr = (1/(muOverRT))*dRhodr
        return [dRhodr, dMdr]


distanceFromCOG = coreRadius+np.linspace(0,numberOfSlices) #number of points
initialDensity = [0,0]       
solution = ode(function,initialDensity,distanceFromCOG).set_integrator('vode', method = 'dopri5')

plt.plot(solution)
solution.set_integrator('dopri5')