Implementing numerical integration in Python

Question

I have this Python code that implements a rectangular numerical integration. It evaluates the (K-1)-dimensional integral for arbitrary integer \$K \geq 1\$

$$\int_{u_K = 0}^{\gamma}\int_{u_{K-1} = 0}^{\gamma-u_K}\cdots\int_{u_2}^{\gamma-u_K-\cdots-u_3}F_U(\gamma-\sum_{k=2}^Ku_k)f_U(u_2)\cdots f_U(u_K)\,du_2\cdots du_K$$

where \$F_U\$ corresponds to the cdf function in the code, and \$f_U\$ to the pdf. I implemented it using recursion as follows:

    #************************** Import necessary libraries******************************************
import numpy as np
import matplotlib.pyplot as plt
import time
#******************************Set the constant scalars and vectors***************************
start_time = time.time()
KU = 3
eta_dB = 3
eta = 10**(eta_dB/10)
ExpanF = 50

tdB = np.arange(-5,11,4)
tVec = 10**(tdB/10)
thVec = (tVec/eta)*(ExpanF-eta*(KU-1))

N = 10000 # For simulation

du = 0.01

#******************************Define functions to be used***************************************
#Define the CDF of U
def CDF(u):
      return 1-1/(u+1)

#Define the PDF of U
def pdf(u):
     return 1/((1+u))**2

def FK(h, k):
    #print(f'h is {h}, and k is {k}')
    if k == 1:
        res = CDF(h)
    else:
        #n_iter = int(h/du)
        res =  0
        u = 0
        while u < h:
            res += FK(h-u, k-1)*pdf(u)*du
            u += du
    return res

#*******************Find the numerical and simulation values of the integral******************
ResultNum = []
ResultSim = []

if (ExpanF-eta*(KU-1)) > 0:
    for t in thVec:
    # Numerical
        ResultNum.append(1-FK(t, KU))
    # Simulation
        count = 0
        for n in range(0, N):
            if np.sum(np.random.exponential(1, (1, KU))/np.random.exponential(1, (1, KU))) <= t:
                count += 1
        ResultSim.append(1-count/N)

The parameter du is the resolution of the numerical integral, and the smaller, the more accurate the numerical approximation is. However, decreasing du comes at the expense of more computational time using the above code. For example I tried to run the above code on my machine Intel i5 @2.5 GH with 6 GB of RAM for du=0.0001, and it took 15 hours and didn't finish, at which case I had to abort it.

Is there anyway the above code can be optimized to run significantly faster? I tried C++, and I got a 30x speed factor using the same code and algorithm, but I was wondering if I can get a similar or better speed up factor in Python.