title: Two_body_problem: scipy.integrate.RK45 gives broadcasting error and scipy.integrate.LSODA never enters the twoBody function
body: <p>I'm working on a trajectory calculator for the Two Body Problem, and I'm attempting to use Scipy's <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html" rel="nofollow noreferrer">RK45</a> or <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.LSODA.html" rel="nofollow noreferrer">LSODA</a> to solve the ODE and return the trajectory. (Please suggest another method if you think there's a better/easier way to do this)</p>
  
  <p>I'm using the Spyder IDE with Anaconda. Scipy version 1.1.0</p>
  
  <p>THE PROBLEMS:</p>
  
  <p>RK45: 
  When using <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html" rel="nofollow noreferrer">RK45</a>, the first step seems to work. When stepping through the code in the debugger, twoBody() is entered, and works exactly as expected the first run through. However, after the first <code>return ydot</code>, things start to go wrong. With a breakpoint at the line <code>ydot[0] = y[3]</code>, we start to see the problem. The array <code>y</code> (which I expected to be 6x1) is now a 6x6 array. When attempting to evaluate this line, numpy returns the error 
  <code>ValueError: could not broadcast input array from shape (6) into shape (1)</code>. Is there an error in my code that would cause <code>y</code> to go from 6x1 to 6x6? Below is the array y right before the broadcasting error is returned.</p>
  
  <pre><code>y = 
   -5.61494e+06   -2.01406e+06    2.47104e+06     -683.979    571.469  1236.76
  -5.61492e+06    -2.01404e+06    2.47106e+06     -663.568    591.88   1257.17
  -5.6149e+06     -2.01403e+06    2.47107e+06     -652.751    602.697  1267.99
  -5.61492e+06    -2.01405e+06    2.47105e+06     -672.901    582.547  1247.84
  -5.61492e+06    -2.01405e+06    2.47105e+06     -672.988    582.46   1247.75
  -5.61492e+06    -2.01405e+06    2.47105e+06     -673.096    582.352  1247.64
  </code></pre>
  
  <p>Could my initial condition <code>Y0</code> be causing it to reach a step too small, and therefore error out?</p>
  
  <p>LSODA:
  I also tried to use the <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.LSODA.html" rel="nofollow noreferrer">LSODA</a> solver. However, it never even enters the <code>twoBody()</code> function! A breakpoint inside at the top of <code>twoBody()</code> is never reached, and the program returns the runtime. I have no idea what's going on here. Guessing I set it up incorrectly.</p>
  
  <p>EDIT:
  The same occurs when using Scipy's <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp" rel="nofollow noreferrer">solve_ivp</a>. All the other methods of integration return the broadcast error.</p>
  
  <pre><code>import numpy as np
  import scipy.integrate as ode
  from time import time
  startTime = time()
  
  def twoBody(t, y):
      """
      Two Body function returns the derivative of the state space variables.
  INPUTS: 
      --- t ---
          A scalar time value. 
  
      --- y ---
          A 6x1 array of the state space of a particle in 3D space  
  OUTPUTS:
      --- ydot ---
          The derivative of y for the two-body problem
  
  """
      mu = 3.986004418 * 10**14
  
      r = np.sqrt(y[0]**2 + y[1]**2 + y[2]**2)
  
      ydot    = np.empty((6,1))
      ydot[:] = np.nan
  
      ydot[0] = y[3]             
      ydot[1] = y[4]             
      ydot[2] = y[5]             
      ydot[3] = (-mu/(r**3))*y[0]
      ydot[4] = (-mu/(r**3))*y[1]
      ydot[5] = (-mu/(r**3))*y[2]
  
      return ydot
  
  
  # In m and m/s
  # first three are the (x, y, z) position
  # second three are the velocities in those same directions respectively
  Y0 = np.array([-5614924.5443320004,
                 -2014046.755686,
                 2471050.0114869997,
                 -673.03650300000004,
                 582.41158099999996,
                 1247.7034980000001])
  
  solution = ode.LSODA(twoBody, t0 = 0.0, y0 = Y0, t_bound = 351.0)
  #solution = ode.RK45(twoBody, t0 = 0.0, y0 = Y0, t_bound = 351.0)
  
  runTime = round(time() - startTime,6)
  print('Program runtime was {} s'.format(runTime))
  </code></pre>