R2 to R2 data fit?

I have some 4 dimensional data. The second two columns (Z and W) are functions of the first two (X and Y). I would like to get a matrix such that

              Z = |a_0    a_1 | X
              W   |a_2    a_3 | Y

is there a good way to do this? I can fit it with a polynomial but then I have a function R1 -> R3 instead of R2 -> R2.