#This sol.py is a script for solving second order ODE's in python. 
from scipy.optimize import minimize
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import draw, show
import decimal
#We want to solve x''+mu*x'+c*x=0. But this method only works
#for systems of first order differential equations. Hence 
# we define 
#x' =: y
#and then 
#y' = -mu*y - c*x 
#The initial values are given by 
#x[0]=-2,x'[0]=y[0]=1

def fun(z,t,mu,c): 
     x,y=z
     return np.array([y, -mu*y - c*x])#[y,y']

mu, c = 4, 4
tmax, dt = 5, 0.01#time increments
t=np.linspace(0, tmax, num=np.round(tmax/dt)+1) 
for mu in [float(j) / 5 for j in range(20, 26)]:
     print 'Solving for mu = {}'.format(mu)
     sol = odeint(fun, [-2, 1], t, args=(mu, c ))[..., 0]
     plt.plot(t, sol, label='mu = {}'.format(mu))
     plt.legend()
show()