import numpy as np
from scipy.integrate import quad,odeint
import matplotlib.pyplot as plt

def J(x) :
    def g(t) :
        return np.cos(x*np.sin(t))
    return quad(g,0,np.pi/2)[0]

def S(x) :
    s = 0
    fact = 1
    for k in range(0,50) :
        s += (-1)**(k)*x**(2*k)/(2**(2*k+1)*fact**2)
        fact *= k+1
    return s*np.pi
    
def f(x,t) :
    a,b = x
    return np.array([b,-b/t-a])
    

X = np.arange(-10,10,0.01)
Y = [J(x) for x in X]
Z = [S(x) for x in X]
plt.plot(X,Y)
plt.plot(X,Z)
X1 = np.arange(0.01,10,0.01)
T = odeint(f,np.array([np.pi/2,0]),X1)
plt.plot(X1,T[:,0])
plt.show()