This seems like a fairly common technique but I'm fairly new to programming and don't quite know what the proper terms to find such an algorithm would be. In my case, it seems like a list of all possible values would approach almost a billion, and I feel like there's a quicker approach to generate the list and match it to a desired value. Essentially, I have a 2D function that I am integrating and I have constants in the function that I want to find all possible values for. Here's the basic integration code: Code (Text): import numpy as np from scipy import integrate i = 0 while i < 1: #run once i = i + 1 sigma0 = 4 r0 = 0 theta0 = 0 def G(r,theta): return (np.e**(-((r**2 + r0**2 - 2*r*r0*(np.cos(theta)*np.cos(theta0) + \ np.sin(theta)*np.sin(theta0)))/(2*sigma0**2))))*r R1I = integrate.nquad(G, [[0,4],[0,2*np.pi]]) R2I = integrate.nquad(G, [[4,7],[0,0.5*np.pi]]) R3I = integrate.nquad(G, [[4,7],[0.5*np.pi,np.pi]]) R4I = integrate.nquad(G, [[4,7],[np.pi,1.5*np.pi]]) R5I = integrate.nquad(G, [[4,7],[1.5*np.pi,2*np.pi]]) In the case above, I specified values for sigma0, r0, and theta0, but I don't want to keep these constant (i.e. I don't want to specify them). What I want to do is set constraints and then find all possible values these can 3 variables can take on. Essentially, I want to specify a range for R1I, R2I, R3I, R4I, and R5I, and then find all possible combination of sigma0, r0, and theta0 in the function that would allow this. Making the range as small as possible and on arbitrary values (e.g. 0.055 < R1I < 0.060, 0.008 < R2I < 0.013, 0.016 < R3I < 0.025, 0.006 < R4I < 0.007, 0.020 < R5I < 0.988) would be optimal. Any advice or recommendations would be greatly appreciated! Thank you!