Transfer function estimation for a MIMO state space system

I have a question regarding the calculation of MIMO transfer functions with the System Identification Toolbox. I would like to calculate the transfer function for a state-space system by using the tfest command. However, I am not sure how to create a system in order to calculate a MIMO System instead of a SIMO. When you run the code and use the “display on” option the toolbox shows:

"Plant Identification Progress"

Transfer Function Identification
Estimation data: Time domain data data
Data has 2 outputs, 2 inputs and 1 408 samples.
Number of poles: [2 2; 2 2],
Number of zeros [1 0; 0 1]

However, the zeros are "swapped". The number of zeros should be as follows: Number of zeros:

 [0 1; 0 1]

In the following code I put on a sample system. I calculate an init system and a system with pole / zero-settings. In my opinion the behavior of "tfSys" estimation is wrong. The representation of the second order system is not right.

Can anyone explain to me, how to create a MIMO state-space-system to work correctly with the tfest command?

Thank you so much!

    % constate values
    a_const = 1; b_const = 1; c_const = 1;

    % systemoutput - simulation
    load t.mat; load x.mat; load sol.mat

    % resampling for the samplingrate
    tnew = 0:.01:t(end);
    y = deval(sol,tnew);

    t = tnew'; x = y';

    %%State-Space-System
    % systemmatrix
    A = [ -1 ,   1;
              1  , 2];

    % inputmatrix
    B = [1, 0;
            0, a_const];  

    % outputmatrix
    C = eye(2);
    D = 0;

    % state-space-system-identification
    init_sys = idss(A, B, C, zeros(size(A)), zeros(size(A)), zeros(size(A,1),1),0);

    % Impuls-Input 2x
    Input_exp = [b_const*ones(size(x,1),1),c_const*ones(size(x,1),1)];

   % simulations-Output 2x
    Output_exp = x;

    % sampling time
    Ts = t(2) - t(1);

    % Timeseries daten
    data = iddata(Output_exp,Input_exp,Ts);

    % Transferfunktions-Modell
    init_sys_TF = idtf(init_sys);

    %%Transferfunction Estimate
    % Transferfunktions-INIT SYSTEM
    tfSys = tfest(data,init_sys_TF,opt,'Ts',Ts); %->wrong 2nd order model

    % Transferfunctions-INIT ZERO/POL
    tfSysPolZero = tfest(data,[2 2;2 2],[0 1; 0 1],opt,'Ts',Ts);

    %%ANALYSE
    figure
    compare(data,tfSys,tfSysPolZero)