numerical ODE solving in python
Question: I am new to Python so at this moment in time, I can only very basic problems.

How do I numerically solve an ODE in Python?

Consider

![equation to solve](http://i.stack.imgur.com/ACR45.gif)

    
    
    \ddot{u}(\phi) = -u + \sqrt{u}
    

with the following conditions

    
    
    u(0) = 1.49907
    

and

    
    
    \dot{u}(0) = 0
    

with the constraint

    
    
    0 <= \phi <= 7\pi.
    

Then finally I want to produce a parametric plot where the x and y coordinates
are generated as a function of u.

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 away to run the twice all at
once.

Here is the error:

    
    
    odepack.error: The function and its Jacobian must be callable functions
    

which the code below generates. The line in question is the sol = odeint.

    
    
    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()
    

**Edit**

I am trying to construct the plot on page 9

[Classical Mechanics
Taylor](http://www.colorado.edu/physics/phys3210/phys3210_fa05/Homework/phys3210_fa05_homework10_solutions.pdf)

Here is the plot with Mathematica

![mathematica plot](http://i.stack.imgur.com/9X6Va.png)

    
    
    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 -> {{-2, 2}, {-2.5, 2.5}}]
    


Answer: 
    import scipy.integrate as integrate
    import matplotlib.pyplot as plt
    import numpy as np
    
    pi = np.pi
    sqrt = np.sqrt
    cos = np.cos
    sin = np.sin
    
    def deriv_z(z, phi):
        u, udot = z
        return [udot, -u + sqrt(u)]
    
    phi = np.linspace(0, 7.0*pi, 2000)
    zinit = [1.49907, 0]
    z = integrate.odeint(deriv_z, zinit, phi)
    u, udot = z.T
    # plt.plot(phi, u)
    fig, ax = plt.subplots()
    ax.plot(1/u*cos(phi), 1/u*sin(phi))
    ax.set_aspect('equal')
    plt.grid(True)
    plt.show()
    

![enter image description here](http://i.stack.imgur.com/lLkco.png)

