I have two fields, `A`

(complex scalar field) and `B`

(real vector field), which are sampled on staggered grids. This means that (in, say, three dimensions) `A`

is an complex array with N^3 elements, while `B`

is a real array with 3*(N-1)^3 elements. I also have corresponding equations of motion with `dA = fA(A,B)`

and `dB = fB(A,B)`

. I’m rolling my own integrator, but I would prefer to switch to DifferentialEquations.jl.

I can see a trivial, even if inconvenient, solution to the difference of array shapes by flattening and stacking both into a vector and then `view`

-ing `reshape`

-ing them into correct form inside the problem definition.

My bigger concern is the complex/real split between the arrays. This type difference has physical relevance and is enforced through the code base, so I would prefer to not cast the real field into a complex-valued field just for sake of this numerical step. There is also space, where in practical situations, N will be large (easily N=10^4, could increase magnitude or two) and I’m not too keen on almost doubling the memory footprint just to store billions of zeros.

What is the most graceful way of solving a system of equations like this?