How can I optimize a 4D integral computation in quantum computing using TensorFlow Quantum and Numpy?

Question

Hello I am new to programming I want to solve this 4 dimensional integral

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fftn, fftshift, ifftn, ifftshift
from numpy import cos, sin, pi

w_0 = 388 * 1e-6 *1e6
Q =28.74 * np.pi/180
L = 2  * 1e-3 *1e6
z = 35 * 1e-3 *1e6
c = 3 * 1e8 *1e6
d = 107.8 * 1e-2 *1e6
theta = Q * 180/np.pi

lambda_p = 405 * 1e-9 *1e6
lambda_s = 810 * 1e-9 *1e6
lambda_i = 810 * 1e-9 *1e6

n_po = np.sqrt(2.7405 + 0.0184/(lambda_p**2 - 0.0179)- 0.0155 * lambda_p**2)

n_pe = np.sqrt(2.3730 + 0.0128/(lambda_p**2 - 0.0179)- 0.0044 * lambda_p**2)

n_so = np.sqrt(2.7405 + 0.0184/(lambda_s**2 - 0.0179)- 0.0155 * lambda_s**2)

n_se = np.sqrt(2.3730 + 0.0128/(lambda_s**2 - 0.0179)- 0.0044 * lambda_s**2)

n_io = np.sqrt(2.7405 + 0.0184/(lambda_i**2 - 0.0179)- 0.0155 * lambda_i**2)

n_ie = np.sqrt(2.3730 + 0.0128/(lambda_i**2 - 0.0179)- 0.0044 * lambda_i**2)

k_s = 2 * np.pi/lambda_s
k_p = 2 * np.pi/lambda_p
k_i = 2 * np.pi/lambda_i
omega_p =  k_p * c
omega_s =  k_s * c
omega_i =  k_i * c

alpha_p = (n_po**2 - n_pe**2) * np.sin(Q) * np.cos(Q)/(n_po**2 * np.sin(Q)**2 + n_pe**2 * np.cos(Q)**2)

beta_p = n_po * n_pe/(n_po**2 * np.sin(Q)**2 + n_pe**2 * np.cos(Q)**2)

gamma_p = n_po/np.sqrt(n_po**2 * np.sin(Q)**2 + n_pe**2 * np.cos(Q)**2)

eta_p = n_po * n_pe/np.sqrt(n_po**2 * np.sin(Q)**2 + n_pe**2 * np.cos(Q)**2)


def f(rho_s,rho_i,phi_s,phi_i,l):
    
    q_sx = rho_s * cos(phi_s)
    q_sy = rho_s * sin(phi_s)
    q_ix = rho_i * cos(phi_i)
    q_iy = rho_i * sin(phi_i)
    q_px = q_sx + q_ix
    q_py = q_sy + q_iy
    
    k_sz =  n_so * omega_s/c - (c/(2 * n_so * omega_s)) * (q_sx**2 + q_sy**2 )
    k_pz =   -alpha_p * q_px + eta_p * omega_p/c 
    k_iz =  n_io * omega_i/c - (c/(2 * n_io * omega_i)) * (q_ix**2 + q_iy**2 )
    Delta_k_z = k_sz + k_iz - k_pz
    
    phi = L * np.sinc(Delta_k_z * L/2)* np.exp(1j * Delta_k_z * L/2)
    V = (w_0/np.sqrt(2 * pi)) * np.exp(-(q_px**2 + q_py**2) * w_0**2/4)            
    return  V * phi * np.exp(1j*l*(phi_s-phi_i))

phi_s = np.linspace(-pi,pi,30)
phi_i = np.linspace(-pi,pi,30)
rho_s = np.linspace(0,0.01,30)
rho_i = np.linspace(0,0.01,30)
l = np.linspace(-50,50,100)

phi_S,phi_I,rho_S,rho_I = np.meshgrid(phi_s,phi_i,rho_s,rho_i)

def summed(rho_s,rho_i,l):
    return np.sum(f(rho_s,rho_i,phi_S,phi_I,l))

def Int(l):
    return np.sum(np.abs(summed(rho_S,rho_I,l))**2*rho_S*rho_I)

result = np.zeros(len(l))
for i, o in enumerate(l):
    result[i] = Int(o)
    
plt.plot(l,result)
plt.show()

I tried solving with nquad also but this 4D integral is taking too much time there. Help me to optimize this code. THanks.I tried solving with nquad also but this 4D integral is taking too much time there. Help me to optimize this code. THanks.I tried solving with nquad also but this 4D integral is taking too much time there. Help me to optimize this code. THanks