mean = [0,0]
low = np.array([-np.inf, -np.inf])
upp = scipy.stats.norm.ppf([0.03, 0.04])
rho = Symbol('rho')
rho = solve((mvn.mvnun(low, upp, mean, [[1,rho],[rho,1]]))[0] - 0.03)
print(rho)