There is no multi-output recursive estimator available today in System Identification Toolbox; this is a limitation. What you need is a true multi-output estimator (e.g., armax, or state-space form), which doesn't exist right now.
Workaround: try some post-processing manipulations like converting to state-space (using the Model Type Converter block), concatenating and possibly balanced reduction (balred).
