Mathematical Graphics

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spi

mu = 1.

def f(u,t):
    dudt = np.array([0.,0.,0.,0.])
    D1 = ((u[0] + 1)**2 + u[2]**2)**(3/2)
    D2 = ((u[0] - 1)**2 + u[2]**2)**(3/2)
    dudt[0] = u[1]
    dudt[1] = -(u[0] + 1)/D1 - mu*(u[0] - 1)/D2
    dudt[2] = u[3]
    dudt[3] = -u[2]/D1 - mu*u[2]/D2
    return dudt

u0 = [0.,1.5,1.,0.]
t = np.linspace(0,200,2000)
u = spi.odeint(f,u0,t)
plt.plot(u[:,0],u[:,2],-1,0,'r*',1,0,'r*'), plt.axis('off')
plt.savefig('euler3body.png',dpi=300)
plt.show()

f = lambda u,t: np.array([1 - u[1]**2,u[0] - u[1]])
t = np.linspace(0,-1,100)
for _ in range(10):
    u0 = [-3,2*np.random.rand() - 3]
    u = spi.odeint(f,u0,t)
    plt.plot(u[:,0],u[:,1],c='b',alpha=0.2)
plt.axis([-3,3,-3,3]), plt.grid(True)
#plt.savefig('phaseportrait.png',dpi=300)
plt.show()