State-Space-Model from more dimensional Transfer-Function
In Matlab, you could just apply ss() to your matrix of Transfer-functions and that should do it.

For example,

G = [tf(rand(1,1),rand(1,2)), tf(rand(1,1),rand(1,2)); tf(rand(1,1),rand(1,2)), tf(rand(1,1),rand(1,2))];

ss(G)


and you get

G

Transfer function 'G' from input 'u1' to output ...

           0.1182
 y1:  -----------------
      0.1398 s + 0.1585

           0.2366
 y2:  -----------------
      0.4479 s + 0.7206

Transfer function 'G' from input 'u2' to output ...

           0.8437
 y1:  ----------------
      0.446 s + 0.4288

           0.8192
 y2:  -----------------
      0.7782 s + 0.0429


ss(G)

ans.a =
               x1          x2          x3          x4
   x1   3.781e-17   5.228e-17  -6.661e-16     0.08869
   x2   6.074e-16   3.508e-16        1.09   9.021e-17
   x3           0          -1      -2.095   9.714e-17
   x4          -1  -5.551e-16           0      -1.664

ans.b =
             u1        u2
   x1  -0.01015   -0.5903
   x2   -0.5066    -1.337
   x3     0.527      1.18
   x4    0.1841    0.3669

ans.c =
          x1     x2     x3     x4
   y1      0      0  1.604      0
   y2      0      0      0  2.869

ans.d =
       u1  u2
   y1   0   0
   y2   0   0

Continuous-time model.