Your question is a bit vague, but let me try the following statement, which might be the kind of answer you are looking for: If $M$ is a manifold and $S\to M$ is a vector bundle over $M$ endowed with a connection $\nabla$, then a 'total differential equation' for sections $s$ of $S$ is an equation of the form $\nabla_X s = 0$ for all vector fields $X$. Then the Frobenius integrability conditions for this system are simply $\Omega^\nabla=0$, where $\Omega^\nabla$ is the curvature of $\nabla$, i.e., it is the $2$-form on $M$ with values in $\mathrm{End}(S)$ that is defined by $$ \Omega^\nabla(X,Y)s = \nabla_X\nabla_Ys - \nabla_Y\nabla_Xs - \nabla_{[X,Y]}s $$ for all vector fields $X$ and $Y$ on $M$. So, yes, the Frobenius integrability conditions are 'covariant', i.e., they can be formulated independently of coordinates. Now, you may be thinking, "But suppose I start with some linear total differential system $$ \frac{\partial u^a}{\partial x^i} = \Gamma^a_{bi}(x) u^b \tag{1} $$ in local coordinates. Does this apply?" The answer is 'yes': You regard the $u^a$ as the components of a section $s = (u^a)$ of a trivialized vector bundle $S$ and the matrices $\Gamma_i = \bigl(\Gamma^a_{bi}(x)\bigr)$ as sections of $\mathrm{End}(S)$ and you define a connection $\nabla$ on this trivial bundle by the rule $$ \nabla_{\frac{\partial}{\partial x^i}}s = \frac{\partial s}{\partial x^i} - \Gamma_i\ s. $$ Then the curvature of $\nabla$ vanishes if and only if the Frobenius conditions are satisfied for the system (1). If you have an inhomogeneous linear system $$ \frac{\partial u^a}{\partial x^i} = \Gamma^a_{bi}(x) u^b + \Lambda^a(x), \tag{2} $$ you can convert it to a homogeneous system by replacing $\Lambda^a(x)$ by $\Lambda^a(x)u^0$ in the above equation and adjoining the equations $$ \frac{\partial u^0}{\partial x^i} = 0, $$ making it a homogeneous linear system in one more variable. Again, the Frobenius conditions are equivalent to the vanishing of the curvature of the connection associated to the augmented system. If you go whole hog to a fully nonlinear total differential system, say, $$ \frac{\partial u^a}{\partial x^i} = \Gamma^a_{i}(x,u), \tag{3} $$ then you can't (usually) reformulate this as the vanishing of the covariant derivatives of a section of something. However, there is still a 'covariant' coordinate-free interpretation of the Frobenius condition in terms of a tensor, but the tensor is defined on the total space $S$ (which has the $x$'s and $u$'s as local coordinates). The point is then, that, on $S$, you have a plane field $D$ (sometimes called a 'distribution') defined by the Pfaffian equations $$ du^a - \Gamma^a_{i}(x,u)\ dx^i = 0 $$ (summation on $i$ intended). Associated to this plane field $D\subset TS$, there is a natural skew-symmetric linear operator $$ \Phi: D\times D \to TS{\bigl/}D $$ i.e., a section of the bundle $\bigl(TS{\bigl/}D\bigr)\otimes \Lambda^2(D^\ast)$ over $S$ that vanishes if and only if the system (3) satisfies the Frobenius conditions. It is this operator that is the 'covariant' interpretation of the Frobenius integrability conditions.