Plotting integral solution with varying parameter

I know there's similar questions around here but none of them get to the root of my problem.

I have an integral which contains a parameter which I want to produce a plot against.

My code is

import numpy as np

import matplotlib.pyplot as plt

from scipy import integrate


def intyfun(x, a):
    return np.exp(-a/4)*np.exp(-x**2*a)*(2*np.pi*x)*(np.sinh(np.pi)/np.cosh(np.pi*x)**2)

Now I'm stuck. I want integrate this function for x over 0 to infinity, and plot the value of it as a varies on the x axis as a parameter. How can I do this?

In mathematica I can do this and the plot looks like this

My mathematica code is

integral[a_?NumericQ] := 
 NIntegrate[
  Exp[-a/4]*Exp[-mu^2*a]*(2*Pi*mu*Sinh[mu*Pi])/(Cosh[mu*Pi]^2), {mu, 
   0, Infinity}, 
  Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0, 
    "MaxErrorIncreases" -> 10000, "SingularityHandler" -> "IMT"}, 
  MaxRecursion -> 100, PrecisionGoal -> 4]

Plot[integral[a], {a, 0.01, 10}, ImageSize -> Large, 
 FrameStyle -> Black,
 BaseStyle -> {FontFamily -> "Latin Modern Roman"}, PlotLabel -> "", 
 PlotStyle -> Black, FrameStyle -> Black,
 BaseStyle -> {FontFamily -> "Latin Modern Roman"}, PlotRange -> All, 
 AxesLabel -> {a, IntegralValue}]

if that helps.

N.B mu=x in my python code.

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


def function(x, a):
    return np.exp(-a/4)*np.exp(-x**2*a)*(2*np.pi*x)*(np.sinh(np.pi)/np.cosh(np.pi*x)**2)


def integrate_function(x_min, x_max, a):
    return integrate.quad(function, x_min, x_max, args=(a,))[0]


# define integration range
x_min = 0
x_max = 1
# define range of a variable
a_points = np.linspace(0, 10, 100)

results = []
for a in a_points:
    value = integrate_function(x_min, x_max, a)
    results.append(value)


plt.plot(a_points, results)
plt.ylabel("Integral value")
plt.xlabel("a variable")
plt.show()

Output: