Matrices form a vector space. Therefore, you can simply integrate them componentwise.
In detail. Let $A:t\mapsto A(t)$ be a function from a real interval $I$ to the space of $m\times n$ real matrices.
Every entry $a_{ij}$ is a real function of a real variable. If all entries are integrable functions, then you can define the integral of the matrix as the matrix of the integrals:
$$\int A(t)\,dt := \left( \int a_{ij}(t)\,dt \right).$$
For your problem:
$$\int \begin{pmatrix}\sin(s)&\cos(\beta s)\\ \cos(s)&\cos(\beta s)\end{pmatrix}ds =
\begin{pmatrix}\int\sin(s)\,ds &\int\cos(\beta s)\,ds\\ \int\cos(s)\,ds&\int\cos(\beta s)\,ds\end{pmatrix} = \begin{pmatrix}-\cos(s)&\frac{1}{\beta}\sin(\beta s)\\ \sin(s)&\frac{1}{\beta}\sin(\beta s)\end{pmatrix}. $$
There is a more sophisticated operation, in case the matrix in question belongs to a Lie algebra: ordered exponentiation. It is to integration as exponentiation is to multiplication, and permits to go from a Lie algebra element (intuitively, a differential transformation) to a group element (a whole transformation).
In this case, you need a $n\times n$ matrix-valued function.
It is explained quite well here:
http://en.wikipedia.org/wiki/Ordered_exponential.
The system is given by
$$\dot{x}=f(x,v)$$
$$\dot{v}=g(x,v).$$
The idea behind the Jacobian is that we use the Taylor expansion of the functions $f$ and $g$ evaluated at the point $x_0$ and $y_0$. The expansions are given by
$$f(x,v)\approx f(x_0,v_0)+\left.\dfrac{\partial f}{\partial x}\right|_{x=x_0,v=v_0}(x-x_0)+\left.\dfrac{\partial f}{\partial v}\right|_{x=x_0,v=v_0}(v-v_0)$$
$$g(x,v)\approx g(x_0,v_0)+\left.\dfrac{\partial g}{\partial x}\right|_{x=x_0,v=v_0}(x-x_0)+\left.\dfrac{\partial g}{\partial v}\right|_{x=x_0,v=v_0}(v-v_0).$$
If $x_0,v_0$ is an equilibrium point we know that $f(x_0,v_0)=g(x_0,v_0)=0$. If we introduce the new state variables $z_1=x-x_0$ and $z_2=v-v_0$ we can rewrite an approximate version of the initial equation as
$$\dot{z}_1=\left.\dfrac{\partial f}{\partial x}\right|_{x=x_0,v=v_0}z_1+\left.\dfrac{\partial f}{\partial v}\right|_{x=x_0,v=v_0}z_2$$
$$\dot{z}_2=\left.\dfrac{\partial g}{\partial x}\right|_{x=x_0,v=v_0}z_1+\left.\dfrac{\partial g}{\partial v}\right|_{x=x_0,v=v_0}z_2.$$
This is a linear system with the system matrix (which is the Jacobian of the nonlinear system) evaluated at $x_0,v_0$.
$$
J=\begin{bmatrix}
\dfrac{\partial f}{\partial x} & \dfrac{\partial f}{\partial v}\\
\dfrac{\partial g}{\partial x} & \dfrac{\partial g}{\partial v}\\
\end{bmatrix}_{x=x_0,v=v_0}
$$
If you permute the columns or the rows nothing will change algebraically, but the correspondence with the variables $z_1$ and $z_2$ will change.
Best Answer
How did you get that $G_v$?
The rest looks good assuming $\omega$ is a constant.