Indeed, this is exactly the kind of situation I wrote xyzpy for. As a quick demo:

```python import numpy as np import xyzpy as xyz

def some_function(x, y, z): return x * np.random.randn(3, 4) + y / z

Define how to label the function's output

runner_opts = { 'fn': some_function, 'var_names': ['output'], 'var_dims': {'output': ['a', 'b']}, 'var_coords': {'a': [10, 20, 30]}, } runner = xyz.Runner(**runner_opts)

set the parameters we want to explore (combos <-> cartesian product)

combos = { 'x': np.linspace(1, 2, 11), 'y': np.linspace(2, 3, 21), 'z': np.linspace(4, 5, 31), }

run them

runner.run_combos(combos) ```

Should produce: ``` 100%|###################| 7161/7161 [00:00<00:00, 132654.11it/s]

<xarray.Dataset> Dimensions: (a: 3, b: 4, x: 11, y: 21, z: 31) Coordinates: * x (x) float64 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 * y (y) float64 2.0 2.05 2.1 2.15 2.2 2.25 2.3 2.35 2.4 2.45 2.5 ... * z (z) float64 4.0 4.033 4.067 4.1 4.133 4.167 4.2 4.233 4.267 4.3 ... * a (a) int32 10 20 30 Dimensions without coordinates: b Data variables: output (x, y, z, a, b) float64 0.6942 -0.3348 -0.9156 -0.517 -0.834 ... ```

And there are options for merging successive, disjoint sets of data (combos2, combos3, ...) and parallelizing/distributing the work.

There are also multiple ways to define functions inputs/outputs (the easiest of which is just to actually return a xr.Dataset), but do let me know if your use case is beyond them or unclear.