issue_comments

It seems that sympy functions returned by lambdfy may not work on arbitrary dimensional inputs, or might not follow broadcasting rules. The workaround might be something like: ```python def vector_func_wrapper(dx, dy, dz, dt): dx, dy, dz, dt = np.broadcast_arrays(dx, dy, dz, dt) # explicitly broadcast args = [a.ravel() for a in [dx, dy, dz, dt]] # convert everything to a vector return vector_funcN(*args).reshape(dx.shape)

xarray.apply_ufunc(vector_func_wrapper, dx, dy, dz, dt) ```

Broadcast allows a variadic number of arguments, so you can write something like: python gx, dy, gz, gt =xr.broadcast(dx, dy, dz, dt)

I'm not very familiar with sympy's lambdify, but I think something like xarray.apply_ufunc(vector_funcN, dx, dy, dz, dt) should work (where dx, dy, dz and dt are all xarray objects, even without manual broadcasting).

If you make each of your arrays xAxisDomian an xarray.DataArray with appropriate dimensions, you can use xarray.broadcast() to convert all of them to common coordinates. xarray.apply_ufunc may also help.

Is this what you want?

```python import numpy as np import xarray as xr

DomainSpaceTimeSize = 5 # using cartesian 4D SpaceTimeDensity = [100, 5] # 100 divisions in space 5 in time

x_coord = np.linspace(-DomainSpaceTimeSize, +DomainSpaceTimeSize, SpaceTimeDensity[0]) y_coord = np.linspace(-DomainSpaceTimeSize, +DomainSpaceTimeSize, SpaceTimeDensity[0]) z_coord = np.linspace(-DomainSpaceTimeSize, +DomainSpaceTimeSize, SpaceTimeDensity[0]) t_coord = np.linspace(0, +DomainSpaceTimeSize, SpaceTimeDensity[1])

scale_func = lambda x, y, z, t: np.cos(1x+2y+3z-4t)

data = scale_func(*np.meshgrid(x_coord, y_coord, z_coord, t_coord))

da = xr.DataArray(data, dims=['x', 'y', 'z', 't'], coords={'x':x_coord, 'y':y_coord, 'z':z_coord, 't':t_coord})

da.sel(x=0.0, y=1.0, z=0.0, t=3.0, method='nearest') <xarray.DataArray ()> array(-0.8039436986070311) Coordinates: x float64 0.05051 z float64 0.05051 t float64 2.5 y float64 0.9596 ```