title: error: function does not have any attributes rk4
body: <p>I have been trying to run the following code for long and every troubleshooting attempt has led me to another. The code is for solving the TOV equations and present a M-R graph for a neutron star. Please if anyone could help with it, it would be greatly appreciated.</p>
  <p>The code has been taken from the following link: <a href="https://www.ictp-saifr.org/schoolgr/Lecture2Creighton.pdf" rel="nofollow noreferrer">https://www.ictp-saifr.org/schoolgr/Lecture2Creighton.pdf</a></p>
  <p>Just to clarify, I have been using Python 3.8, and the code is 9 years old so that might as well be the problem but i tried running it on a 2.x compiler and faced the same issue. Thanks in advance.</p>
  <pre><code>from scipy.integrate import odeint
  import numpy as np
  import matplotlib.pyplot as plt 
  from scipy.constants import pi, G, c, hbar, m_n
  
  Msun = 1.98892e30
  
  Gamma0 = 5.0/3.0
  K0 = (3.0*pi**2)**(2.0/3.0)*hbar**2/(5.0*m_n**(8.0/3.0))
  Gamma1 = 3.0
  rho1 = 5e17
  P1 = K0*rho1**Gamma0
  K1 = P1/rho1**Gamma1
  
  def eos(rho):
      if rho &lt; rho1: return K0*rho**Gamma0
      else: return K1*rho**Gamma1
  
  def inveos(P):
      if P &lt; P1: return (P/K0)**(1.0/Gamma0)
      else: return (P/K1)**(1.0/Gamma1) 
  
  def tov(y,r):
      P, m = y[0], y[1]
      rho = inveos(P)
      dPdr = -G*(rho + P/c**2)*(m + 4.0*pi*r**3*P/c**2)
      dPdr = dPdr/(r*(r - 2.0*G*m/c**2))
      dmdr = 4.0*pi*r**2*rho
      return np.array([dPdr, dmdr])
  
  def tovsolve(rhoc):
      r = np.arange(10.0, 20000.0, 1)
      m =  np.zeros_like(r)
      P =  np.zeros_like(r)
      m[0] = 4.0*pi*r[0]**3*rhoc
      P[0] = eos(rhoc)
      y = np.array((P[0],m[0]))
      i = 0
      while P[i] &gt; 0.0 and i &lt;len(r)-1:
          dr = r[i+1] - r[i]
          y = odeint.rk4(tov, y, r[i], dr)
          P[i+1] = y[0]
          m[i+1] = y[1]
          i = i+1
      return m[i-1]/Msun, r[i-1]/1000
  
  rhoc = np.logspace(17.5,20)
  M = np.zeros_like(rhoc)
  R = np.zeros_like(rhoc)
  for i in range(len(rhoc)):
      M[i], R[i] = tovsolve(rhoc[i])
  
  plt.plot(R,M)
  plt.xlabel('Radius (km)')
  plt.ylabel('Mass (Solar Masses)')
  plt.grid()
  plt.show() 
  </code></pre>