Как ускорить вычисления построения сложной функции в python?

Question

У меня есть функция f (r, Q), где r и Q кратны 100. Я хочу создать график contour3D для этой функции для всех комбинаций r и Q, (r на оси X, Q по оси Y и f (r, Q) по оси Z.), где 0

import matplotlib.pyplot as plt
import numpy as np
from numpy import sqrt, pi, exp
from scipy.integrate import quad
from mpl_toolkits.mplot3d import Axes3D


mean, std = 295, 250
l, m, p = 7, 30, 15
w, K, c = 7, 100, 5
h = 0.001  # per unit per  day

# defining Cumulative distribution function
def cdf(x):
    cdf_eqn = lambda t: (1 / (std * sqrt(2 * pi))) * exp(-(((t - mean) ** 2) / (2 * std ** 2)))
    cdf = quad(cdf_eqn, -np.inf, x)[0]
    return cdf
# defining Probability density function
def pdf(x):
    return (1 / (std * sqrt(2 * pi))) * exp(-(((x - mean) ** 2) / (2 * std ** 2)))
# getting the equation in place
def G(r, Q):
    return K + c * Q \
           + w * (quad(cdf, 0, Q)[0] + quad(lambda x: cdf(r + Q - x) * cdf(x), 0, r)[0]) \
           + p * (mean * l - r + quad(cdf, 0, r)[0])
def CL(r, Q):
    return (Q - r + mean * l - quad(cdf, 0, Q)[0]
            - quad(lambda x: cdf(r + Q - x) * cdf(x), 0, r)[0]
            + quad(cdf, 0, r)[0]) / mean
def I(r, Q):
    return h * (Q + r - mean * l - quad(cdf, 0, Q)[0]
                - quad(lambda x: cdf(r + Q - x) * cdf(x), 0, r)[0]
                + quad(cdf, 0, r)[0]) / 2
# pulling it all together
def f(r, Q):
    TC = G(r, Q)/CL(r, Q) + I(r, Q)
    return TC

fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111, projection='3d')

'''need to change the step size to 100 and do,1000 is for 
the faster computation and it may take a while, you can experiment with 
different step size values to make it faster, try 2000, 3000 as step size
to make it computationally faster'''

r = list(range(0,5000,1000))
q = list(range(2000,15000,1000))
R, Q = np.meshgrid(r, q)

zs = np.array([f(r,q) for r,q in zip(np.ravel(R), np.ravel(Q))])
Z = zs.reshape(R.shape)
ax.contour3D(R, Q, Z, 50, cmap='binary')

#ax.plot_surface(R, Q, Z)

ax.set_xlabel('R Label')
ax.set_ylabel('Q Label')
ax.set_zlabel('Z Label')

plt.show()