def METHOD_NAME(n_reps=100, n_extra_axes=1):
    reps = np.arange(n_reps)
    axes = ['sample'] + [f'ax_{i}' for i in range(n_extra_axes)]
    extra_axes_vals = [np.linspace(-1., 1., 10+i) for i in range(n_extra_axes)]
    axes_vals = [reps] + extra_axes_vals
    axes_vals_mesh = np.meshgrid(*axes_vals, indexing='ij')
    outcomes_mesh = np.random.normal(
        loc=0, scale=1, size=axes_vals_mesh[0].shape,
    )
    for axvals in axes_vals_mesh[1:]:
        outcomes_mesh = np.add(outcomes_mesh, axvals)
    dataset = DataDict(
        result=dict(values=outcomes_mesh.flatten(),
                    axes=axes),
    )
    for ax, axvals in zip(axes, axes_vals_mesh):
        dataset[ax] = dict(values=axvals.flatten())
    return dataset