Python - How to quicken stochastic differential equations? Random generation is slow

I am integrating numerically a set of differential equations with noise (stochastic differential equations). The problem is that the noise, which I generate with random normal numbers, is in a for loop, which makes it take forever. This is my code:

def Stochastic_Bloch_Equations(t, v):  #Stochastic differential equation
    return [-0.5*gamma*v[0]-2*(v[3]+x_c)*v[2],\
            -0.5*gamma*v[1]-2*(v[4]+y_c)*v[2],\
            -gamma*v[2]+0.5*gamma+2*v[3]*v[0]+2*v[4]*v[1],\
            -b_x*v[3]+cmth.sqrt(0.5*gamma*(Num_Photons(E, gamma)+Corr_Photons(E, gamma))*b_x**2)*np.random.randn(),\
            -b_y*v[4]+cmth.sqrt(0.5*gamma*(Num_Photons(E, gamma)-Corr_Photons(E, gamma))*b_y**2)*np.random.randn()]

#Parameters
Tmax=100 ##max value for time
Nmax=10000 ##number of steps
dt=Tmax/Nmax  ##increment of time

t=np.linspace(0.0,Tmax,Nmax+1)

g_c=np.linspace(0.1,100.0,20)

sz_c=[]

for gamma_c in g_c:
    b_x=0.5*gamma_c-E
    b_y=0.5*gamma_c+E
    stochastic_sol=solve_ivp(Stochastic_Bloch_Equations, [min(t),max(t)],v_0, t_eval= t) #solver
    sz_c.append(stochastic_sol.y[2,-1]) #take last value for sz, steady state solution
    print(gamma_c)

Notice how for every gamma_c, I integrate numerically the stochastic differential equations at a number of time steps Nmax=10001, which means that it generates 20002 random numbers (I call the random.randn() twice). At the end of the for loop, I would have generated 400040 random numbers. And the whole process is dreadfully slow.

I thought I might generate a vector R with all the necessary random numbers, which is infinitely much faster, and then use the elements R[i] in my function. This way all the random numbers are already generated and I don't have to generate them for every loop and every dt in the integration.

Could anybody guide me as how can I adapt this random vector R such that each element R[i] is called in my function be integrated?

Thanks,