title: Numerical ODE solving in Python
body: <p>How do I numerically solve an ODE in Python?</p>
  
  <p>Consider </p>
  
  <p><img src="https://i.stack.imgur.com/ACR45.gif" alt="equation to solve"></p>
  
  <pre><code>\ddot{u}(\phi) = -u + \sqrt{u}
  </code></pre>
  
  <p>with the following conditions </p>
  
  <pre><code>u(0) = 1.49907
  </code></pre>
  
  <p>and</p>
  
  <pre><code>\dot{u}(0) = 0
  </code></pre>
  
  <p>with the constraint </p>
  
  <pre><code>0 &lt;= \phi &lt;= 7\pi.
  </code></pre>
  
  <p>Then finally, I want to produce a parametric plot where the x and y coordinates are generated as a function of u.</p>
  
  <p>The problem is, I need to run odeint twice since this is a second order differential equation.
  I tried having it run again after the first time but it comes back with a Jacobian error.  There must be a way to run it twice all at once.</p>
  
  <p>Here is the error:</p>
  
  <blockquote>
    <p>odepack.error: The function and its Jacobian must be callable functions</p>
  </blockquote>
  
  <p>which the code below generates.  The line in question is the sol = odeint.</p>
  
  <pre><code>import numpy as np
  from scipy.integrate import odeint
  import matplotlib.pyplot as plt
  from numpy import linspace
  
  
  def f(u, t):
      return -u + np.sqrt(u)
  
  
  times = linspace(0.0001, 7 * np.pi, 1000)
  y0 = 1.49907
  yprime0 = 0
  yvals = odeint(f, yprime0, times)
  
  sol = odeint(yvals, y0, times)
  
  x = 1 / sol * np.cos(times)
  y = 1 / sol * np.sin(times)
  
  plot(x,y)
  
  plt.show()
  </code></pre>
  
  <p><strong>Edit</strong></p>
  
  <p>I am trying to construct the plot on page 9</p>
  
  <p><a href="http://www.colorado.edu/physics/phys3210/phys3210_fa05/Homework/phys3210_fa05_homework10_solutions.pdf" rel="nofollow noreferrer">Classical Mechanics Taylor</a></p>
  
  <p>Here is the plot with Mathematica</p>
  
  <p><img src="https://i.stack.imgur.com/9X6Va.png" alt="mathematica plot"></p>
  
  <pre class="lang-none prettyprint-override"><code>In[27]:= sol = 
   NDSolve[{y''[t] == -y[t] + Sqrt[y[t]], y[0] == 1/.66707928, 
     y'[0] == 0}, y, {t, 0, 10*\[Pi]}];
  
  In[28]:= ysol = y[t] /. sol[[1]];
  
  In[30]:= ParametricPlot[{1/ysol*Cos[t], 1/ysol*Sin[t]}, {t, 0, 
    7 \[Pi]}, PlotRange -&gt; {{-2, 2}, {-2.5, 2.5}}]
  </code></pre>