I guess at this point it would be good if you could go into more detail about your exact model equations? You started from:
I made some small changes from your original equation, because I don’t make a distinction between X and U. They would just be concatenated together into what I call X, while matrices you called A and B would be block concatenated what I called A. I made some simple choices for A, C, Q, but they’re not necessarily appropriate at all. Actually, you can accomplish what I did much easier just using X_t \sim \text{GaussianRandomWalk}(\cdot).
So what kind of dynamics are you modeling? ARIMA? Seasonal cycles? Random walk with drift? Exponential smoothing?
For multivariate, the original model is multivariate. You shouldn’t be using a 3d tensor as the observed data, you just want to add more columns to X, y, or both.