matrix differential equation support

[ArnoStrouwen]

The documentation says CuArray u0 are supported.
But only one dimensional arrays are supported because of lines like:

DiffEqGPU.jl/src/DiffEqGPU.jl

Line 345 in 8489a3e

u0 = reduce(hcat, Array(probs[i].u0) for i in 1:length(I))

To generalize this to N dimensions the u0 should concatenate across an additional dimension.
And the kernels should then index over this last dimension, instead of assuming there are only 2:

DiffEqGPU.jl/src/DiffEqGPU.jl

Line 24 in 8489a3e

@views @inbounds f(du[:, i], u[:, i], p[:, i], t)