Let $X'=AX + b(t)$ be a system of differential equations where: $X(t) = (x_1(t),...,x_n(t))$, $A\in\mathcal M_n({\mathbb{R}})$ and $b(t) = (b_1(t),..,b_n(t))$.
How do I solve this system using the variations of constants method?
I know that if $M(t)=(x_1(t)|....|x_n(t))$ a matrix of fundamental solutions and
$$X(t)=c_1(t)x_1(t)+...+c_n(t)x_n(t)=M(t)C(t)$$
where $C(t) = (c_1(t)|..|c_n(t))$
And that $C(t)$ verifies $C'(t)=M^{-1}(t)b(t)$
If I have an example like:
$$\left\{\begin{array}{l}x'=x+y-\cos t\\y'=-2x-y+\sin t+\cos t\end{array}\right.$$
How do I actually solve it?
You were almost done, as you can simply integrate the function $M^{-1}(t)b(t)$ over $t$ to obtain the function $C(t)$, up to constants of integration.
You write your example in the form \begin{equation} \boldsymbol{\dot{x}} = \left( \begin{array}{c} \dot{x}\\ \dot{y} \end{array} \right) = \left( \begin{array}{cc} 1 & 1\\ -2 & -1 \end{array} \right) \left( \begin{array}{c} x\\ y \end{array} \right) + \left( \begin{array}{c} -\cos(t)\\ \sin(t)+\cos(t) \end{array} \right) = \boldsymbol{\underline{A}} \boldsymbol{x} + \boldsymbol{b}(t), \end{equation} which has the general solution \begin{equation} \boldsymbol{x}(t) = \boldsymbol{\underline{M}}(t) \int \boldsymbol{\underline{M}}^{-1}(t) \boldsymbol{b}(t) \, \mathrm{d}t, \end{equation} where $\boldsymbol{\underline{M}}(t) := e^{t \boldsymbol{\underline{A}}}$ is a matrix exponential (with inverse $\boldsymbol{\underline{M}}^{-1}(t) = e^{-t \boldsymbol{\underline{A}}}$). The constants of integration are contained in the indefinite integral. We obtain \begin{equation} \boldsymbol{\underline{M}}(t) = e^{t\boldsymbol{\underline{A}}} = \left( \begin{array}{cc} \cos(t) + \sin(t) & \sin(t)\\ -2\sin(t) & \cos(t) - \sin(t) \end{array} \right), \quad \int \boldsymbol{\underline{M}}^{-1}(t) \boldsymbol{b}(t) \, \mathrm{d}t = \left( \begin{array}{c} C_1 - t\\ C_2 + t \end{array} \right), \end{equation} $C_1, C_2 \in \mathbb{R}$, and finally \begin{equation} \boldsymbol{x}(t) = \left( \begin{array}{c} x(t)\\ y(t) \end{array} \right) = \left( \begin{array}{c} -t \cos(t) + C_1 (\cos(t) + \sin(t)) + C_2 \sin(t)\\ t \cos(t) + t \sin(t) - 2 C_1 \sin(t) + C_2 (\cos(t)-\sin(t)) \end{array} \right). \end{equation}