The idea is using a "covariant" formalism in classical mechanics. I do not know how accommodate your equations in particular, but a good notion of covariant derivative that physically encompasses the "inertial forces" (as in GR) is at disposal also in classical physics.
Consider the classical spacetime $M$. There there is a privileged class of reference frames called inertial such that changing (Cartesian) coordinates among them
$$t=t'+c\:, \quad x^{i} = R^i\:_j x^{'j} + v^it' + b^i\:, \quad i,j= 1,2,3 \tag{1}$$
for $R\in O(3)$, $c, v^i, b^i\in \mathbb R$
Instead, general transformations of coordinates (always Cartesian) at rest of generic reference frames are all of the form (in view of the fact that lengths are absolute at every given time)
$$t=t'+c\:, \quad x^{i} = R^i\:_j(t') x^{'j} + b^i(t')\:, \quad i,j= 1,2,3 \tag{2}$$
for arbitrary (smooth) functions $\mathbb R \ni t \mapsto R(t) \in O(3)$ and
$\mathbb R \ni t \mapsto b^i(t) \in O(3)$.
In inertial reference frames inertial motions are described by geodesics
$$t= v^0s+ c^0\:, \quad x^i(s) = sv^i + c^i\:,\quad v^0 \neq 0\:.$$
Indeed, they these are solutions $\gamma : \mathbb R \to M$ of the geodesic equation
$$\frac{d^2x^a}{ds^2} + \Gamma^{a}_{bc}\frac{dx^b}{ds}\frac{dx^c}{ds}=0\:, \quad a,b,c =0,1,2,3$$
where $$\Gamma^{a}_{bc}=0\:,\tag{3}$$
with the requirement that $v^0\neq 0$, i.e. $\langle dT, \dot{\gamma}\rangle \neq 0 $, where $T: M \to \mathbb R$ is the absolute classical time, defined up to an additive constant and always coinciding with the coordinate $t$ and $t'$ (up to an additive constant).
All that means that the spacetime of classical physics, as soon as we switch on dynamics according to Newtonian picture, a non-metric, symmetric affine connection arises determining inertial motions.
This connection is trivial just in the Cartesian coordinates at rest with every inertial reference frame. Indeed, it is easy to see that the condition (3) is invariant under transformation of coordinates of the form (1).
Instead, passing to a non-inertial reference frame equipped with (Cartesian) coordinates $x^{'a}$ where $x^{'0}=t'$ by means of (2), we have
$$\Gamma^{'a}_{bc} \neq 0$$
and, referring to (2) where the coordinates $t,x^i$ are supposed at rest with an inertial reference frame, since
$$\Gamma^{'a}_{bc} =\frac{\partial x^{'a}}{\partial x^d} \frac{\partial^2 x^{d}}{\partial x^{'b}\partial x^{'c}} + 0\:, $$
we have
$$\Gamma^{'0}_{bc} =0\:, \quad \Gamma^{'i}_{0j}= \Gamma^{'i}_{j0} = (R^{-1})^{i}_{m}\dot{R}^m_j(t)\:,\quad \Gamma^{'i}_{00}= (R^{-1})^{i}_{m}\ddot{R}(t)^{m}_jx^{'j} + (R^{-1})^{i}_{m}\ddot{b}^m(t)\:, \quad \Gamma^{'i}_{jk}=0\:.$$
In the non-inertial reference frame the geodesic equation (inertial motion) is the complete geodesic equation, where I have introduced the inertial mass $m$ of the point of matter considered,
$$m\frac{d^2x^{'a}}{ds^2}= -m\Gamma^{'a}_{bc}\frac{dx^{'b}}{ds}\frac{dx^{'c}}{ds}\quad a,b,c =0,1,2,3$$
and the terms $-\left(\frac{dx^{'0}}{ds} \right)^{-2} m\Gamma^{'i}_{bc}\frac{dx^{'b}}{ds}\frac{dx^{'c}}{ds}$ are nothing but the various inertial forces written into a more abstract formalism.
Equipped with this affine connection (and taking advantage of the absolute spatial metric $\delta_{ij}$), every differential equation can be written in every coordinate system using the standard procedure, replacing every derivative $\partial_a$ with the corresponding covariant derivative
$\nabla_a$. Inertial forces are automatically embodied in the formalism.
The technically difficult (but not impossible obviously) point is re-writing the coefficients $\Gamma$ in terms of $\vec{\omega}$ its derivative and the acceleration of the origin of the inertial reference frame with respect to the inertial reference frame. By comparing with the classical formula of inertial forces, if $\vec{\omega} = \omega^i {\bf e}'_i$ is the angular velocity of the non-inertial reference frame with respect to the inertial one, $\vec{x'} = x^{'i}{\bf e}'_i$ is the position vector in the non-inertial reference frame and $\vec{a}_{O'} = a^{'i}{\bf e}'_i$ is the acceleration of the origin of axes of the non-inertial reference frame with respect to the inertial one, and everything is represented in the axes ${\bf e}'_1, {\bf e}'_2, {\bf e}'_3$ of the non-inertial frame, it should be (using the Euclidean metric $\delta_{ij}$ to raise and lower indices)
$$\Gamma^{'i}_{0j} = \Gamma^{'i}_{j0} = \omega^k \epsilon^i\:_{jk}\quad \mbox{Coriolis' force}$$
and
$$\Gamma^{i}_{00} = a_{O'}^i + \epsilon^i\:_{jk}\dot{\omega}^j x'^{k} + \omega^i x^j\omega_j - x^i \omega^j \omega_j \quad \mbox{drag force}\:.$$
I am not completely sure on this result as is the first time I try to get it.
