How do I get the transfer functions of each channel for a control system?
This is a good question. If you have a MIMO system like you have described:
$$G(s)=\begin{pmatrix}
G_{11}(s) &amp; G_{12}(s)\\ 
 0&amp;G_{22}(s) 
\end{pmatrix}$$
And if you want to evaluate the channels separately it gets a bit tricky. For the channel \$y_1(s)\$ you certainly cannot do:
$$y_{1}(s)=\frac{G_{11}(s)+G_{12}(s)}{1+G_{11}(s)+G_{12}(s)}$$
This would be true only for the inputs \$u_1 = u_2\$. Therefore you will have to evaluate the influence of the each input separately, what due to the linearity of the model should be general enough. 
$$ y_{1}(s)=\frac{G_{11}(s)}{1+G_{11}(s)}u_1(s)+\frac{G_{12}(s)}{1+G_{12}(s)}u_2(s) $$

But it depends of course what kind of analysis do you want to perform.
In matalb, for MIMO system you can use matrix representation directly.

G11 = tf(...)
G12 = tf(...)
G22 = tf(...)
G = [G11, G12; 0, G22]


This will give you the two input two output MIMO system. To better see the coupling, it is common to continue in state space representation for MIMO systems and in that case you can just use:

G = ss([G11, G12; 0, G22])


To see the poles of each channel you can check the poles of each row of the matrix:

pole(G(1,:)); % y_1
pole(G(2,:)); % y_2


But pole,bode,step, pzmap and similar functions will always give you superposed solutions for all the transfer functions (G11 and G12).

If you really want to dig deeper in this topic I suggest you to check singular values decomposition of the dynamical systems in Matlab you can calculate them using the sigma function. The output of the sigma function is a frequency dependent worst case system gain for each input, taking in consideration the coupling. This is a type of bode plot of the system and you can use the same evaluation tools as for every other bode plot.
In fact for SISO system singular values plot is exactly bode magnitude plot. So in your case your singular value plot will have two lines (due to two inputs).
With the singular value plot you will be able determine the worst case gains and magnitudes on certain frequencies, but not really the poles and zeros of the system.   

EDIT: Control strategy explanation

First of all this equation is not correct:

F4_cross = (F4_min(1,1)/(1+F4_min(1,1)))+(F4_min(1,2)/(1+F4_min(1,2)));


As explained earlier you cannot simply add your transfer functions. 
In order to study this better you need a block diagram like this one:



Furthermore you can write this system as two separate systems:


And the one which has coupling influence, in this case it is just the disturbance for this system.


When you write it this way you can clearly see that the there are two loops you are interested to influence with your contorller:
$$ G_{cl} = \frac{y_1}{r_1} = \frac{C_{1}G_{11}}{1+ C_{1}G_{11}} $$
This is the control loop in which you need to have zero error, therefore you will need to add the integrator if necessary and so on... 

Then, the second transfer funciton you are interested in is:
$$ G_{dist}(s) = \frac{y_1}{u_2} = \frac{G_{12}}{1+C_{1}G_{11}} $$ 
Now in order for your tracking error to be zero you will need to influence this transfer function to reject the disturbance (to diminish the signal \$u_2\$). Meaning that this transfer function should have stationary gain of 0.
So just study both the transfer function \$G_{cl}\$ and \$G_{dist}\$ and adapt controller \$C1\$ if necessary.