Integral of the secant function

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

sec = lambda x: 1./np.cos(x)

# evaluate the secant integral for -pi/2 < theta < pi/2
n = 100
theta = np.linspace(-np.pi/2+0.01, np.pi/2-0.01, n)

# Numerical answer by quadrature
I = np.zeros(n)
for i, th in enumerate(theta):
    I[i], _ = quad(lambda x: sec(x), 0, th)

# Analytical answer: the inverse Gudermannian function
gdinv = np.log(np.abs(sec(theta) + np.tan(theta)))

plt.plot(theta, I, 'o', label='Numerical answer')
plt.plot(theta, gdinv, label='Analytical answer')
plt.legend(loc=4)
plt.show()