Typically with a multiple input, multiple output (MIMO) system, a control engineer uses a state feedback controller. This style of controller leverages a state-space model of the system and generally takes the form:
$$
\dot{x}=\mbox{A}x+\mbox{B}u \\
y = \mbox{C}x + \mbox{D}u \\
$$
where $x$ is a vector of states, $u$ is a vector of inputs, $y$ is a vector of outputs, and the time derivative of the states, $\dot{x}$, shows how the states evolve over time, as determined by combinations of states $\mbox{A}$ and inputs $\mbox{B}$. Outputs are also determined by an interaction between states and inputs, but the outputs can be any combination, so the output state and input matrices are different - $\mbox{C}$ and $\mbox{D}$.
I won't go into a large amount of detail regarding state feedback controls, but in general, the matrices $\mbox{A} \rightarrow \mbox{D}$ "map" or associate a particular state or input to another state or input. For instance, if you want to model a system of unrelated differential equations, you would get something like:
$$
\dot{x} = \left[ \begin{array}{ccc}
\dot{x}_1 \\
\dot{x}_2 \\
\dot{x}_3 \end{array} \right] = \left[ \begin{array}{ccc}
k_1 & 0 & 0 \\
0 & k_2 & 0 \\
0 & 0 & k_3 \end{array} \right] \left[ \begin{array}{ccc}
x_1 \\
x_2 \\
x_3 \end{array} \right]
$$
which represents:
$$
\dot{x}_1 = k_1 x_1 \\
\dot{x}_2 = k_2 x_2 \\
\dot{x}_3 = k_3 x_3 \\
$$
If you wanted to add input $u_1$ to the equation for $\dot{x}_1$ and input $u_2$ to $\dot{x}_3$, then you could add a $\mbox{B}u$ term:
$$
\dot{x} = \left[ \begin{array}{ccc}
\dot{x}_1 \\
\dot{x}_2 \\
\dot{x}_3 \end{array} \right] = \left[ \begin{array}{ccc}
k_1 & 0 & 0 \\
0 & k_2 & 0 \\
0 & 0 & k_3 \end{array} \right] \left[ \begin{array}{ccc}
x_1 \\
x_2 \\
x_3 \end{array} \right] + \left[ \begin{array}{ccc}
1 & 0 \\
0 & 0 \\
0 & 1 \end{array} \right] \left[ \begin{array}{ccc}
u_1 \\
u_2 \end{array} \right]
$$
If you want to keep this, but you think that state $x_1$ contributes to how $x_2$ changes, you can add that interaction:
$$
\dot{x} = \left[ \begin{array}{ccc}
\dot{x}_1 \\
\dot{x}_2 \\
\dot{x}_3 \end{array} \right] = \left[ \begin{array}{ccc}
k_1 & 0 & 0 \\
\boxed{ k_{x_1 \rightarrow x_2} } & k_2 & 0 \\
0 & 0 & k_3 \end{array} \right] \left[ \begin{array}{ccc}
x_1 \\
x_2 \\
x_3 \end{array} \right] + \left[ \begin{array}{ccc}
1 & 0 \\
0 & 0 \\
0 & 1 \end{array} \right] \left[ \begin{array}{ccc}
u_1 \\
u_2 \end{array} \right]
$$
When you write these out now, you get:
$$
\begin{array} \\
\dot{x}_1 & = & k_1 x_1 + u_1 \\
\dot{x}_2 & = & k_{x_1 \rightarrow x_2}x_1 + k_2 x_2 \\
\dot{x}_3 & = & k_3 x_3 + u_2 \end{array}
$$
You can keep building up complexity as your system requires. Once you have a model, for state feedback controls, you need to make sure that the system is linear, in that the system doesn't have trig functions or one state multiplying itself or another state, and make sure that it is time invariant, in that the matrices $\mbox{A} \rightarrow \mbox{D}$ don't change with time - no function of (t) in them. You may be able to make some simplifications, such as a small angle approximation to help get your $\mbox{A}$ matrix into the LTI form required for the next step.
Now you can "mask" the entire system into the tidy two equations first shown, hiding the entire $\mbox{A}$ matrix with just the letter 'A', etc. With the Laplace transform you can (hand-wave) evaluate the uncontrolled, open-loop dynamics of the system. You do this by finding the poles of the system, which in term indicate system response.
You can also evaluate the system to see if it is controllable, meaning that you can use your inputs to alter all of the states in a unique manner, and to see if it is observable, meaning that you can actually determine what the values of the states are.
If the system is controllable, you can take information about the states, $-\mbox{G}x$, and feed that into the system, using the information you have about the states to drive them to a desired value. Using only the two initial equations for clarity, when you add the control signal to the input you get:
$$
\dot{x} = \mbox{A}x + \mbox{B}(u - \mbox{G}x) \\
y = \mbox{C}x + \mbox{D}u \\
$$
which becomes:
$$
\dot{x} = \mbox{A}x - \mbox{BG}x + \mbox{B}u \\
y = \mbox{C}x + \mbox{D}u \\
$$
which can be rearranged as:
$$
\dot{x} = [\mbox{A}-\mbox{BG}]x + \mbox{B}u \\
y = \mbox{C}x + \mbox{D}u \\
$$
Where before you system response was driven by the $\mbox{A}$ matrix, now it is driven by $\mbox{A-BG}$. You can again evaluate the poles via the Laplace transform, but now you have a gain matrix $\mbox{G}$ you can use to tune the controller, putting the poles wherever you want, which establishes the time response to be whatever you want.
The process continues, with observers setup to compare the actual system output $y$ with the model's predicted output $\hat{y}$. This is where it's important to note that the outputs don't have to be the same combination of states as you use in the state differential equation - where your states might be a current your output might be a voltage ($R\times I$) so you can make a comparison with a measurable signal on your real system.
Like I said, there is a ton of information involved with modeling systems and designing state feedback controllers, I just outlined the general process as I believe this is scope you were looking for with your question.
