Question

Jeanne's calculations give the right answer, i.e., that the solutions depend on two arbitrary functions of 2 variables.

It turns out, though, that, with the right choice of variables, one can reduce the problem to an underdetermined system of $2$ linear equations for a pair of tensors of rank $2$ on the $2$-sphere. (This choice of variables is suggested by examining the characteristic variety and tableau of the EDS.)

First of all, if one solves for the coefficients of the $U^i$ as the OP suggests, one does get $7$ parameters. One way to do this is as follows: $$ \begin{align} U^1 &= -(\sqrt{3}/3)\bigl((v_4{-}(\cosh(\rho)-\sinh(\rho))v_1)\ \theta^1 -(v_2{+}(\cosh(\rho)+\sinh(\rho))v_3)\ \theta^2\bigr)\\ U^2 &= (\cosh(\rho)-\sinh(\rho))v_1\ \theta^1-(\cosh(\rho)+\sinh(\rho))v_3\ \theta^2\\ U^3 &= (\cosh(\rho)-\sinh(\rho))v_5\ \theta^1 + v_6\ \theta^2 + 2\sinh(\rho)v_3\ \theta^3 +v_4\ d\rho\\ U^4 &= v_7\ \theta^1-(\cosh(\rho)+\sinh(\rho))v_5\ \theta^2 + 2\sinh(\rho)v_1\ \theta^3 +v_2\ d\rho\\ U^5 &= (\cosh(\rho)-\sinh(\rho))v_2\ \theta^1-(\cosh(\rho)+\sinh(\rho))v_4\ \theta^2.\\ \end{align} $$ The $v_i$ are now $7$ unknowns, and substituting this into the above $2$-forms gives an involutive differential system for the $v_i$ that is generated by five $2$-forms and has characters $s_1=5$ and $s_2=2$, with $s_i=0$ for $i>2$.

However, examining the characteristic variety and tableau of this system suggests that one should reparametrize the first four $v_i$ as $$ \begin{align} v_1 &= -(2\sinh(\rho){+}\cosh(\rho))\ p_2-\cosh(\rho)\ q_2\\ v_2 &= (2\cosh(\rho)\sinh(\rho)+2\cosh^2(\rho)-1)(p_1{+}q_1)\\ v_3 &= +(2\sinh(\rho)-\cosh(\rho))\ p_1 + \cosh(\rho)\ q_1\\ v_4 &= (2\cosh(\rho)\sinh(\rho)-2\cosh^2(\rho)+1)(p_2{-}q_2).\\ \end{align} $$ for some functions $p_1,p_2,q_1,q_2$. When one does this, one finds that the differential equations imply $$ \begin{align} dp_1 &= (p_0{+}p_3)\ \theta^1 + (p_4{+}p_5)\ \theta^2-3p_2\ \theta^3\\ dp_2 &= (p_4{-}p_5)\ \theta^1 + (p_0{-}p_3)\ \theta^2+3p_1\ \theta^3\\ dq_1 &= (q_0{+}q_3)\ \theta^1 + (q_4{+}q_5)\ \theta^2-1q_2\ \theta^3\\ dq_2 &= (q_4{-}q_5)\ \theta^1 + (q_0{-}q_3)\ \theta^2+1p_1\ \theta^3\\ \end{align} $$ for some functions $p_0,p_3,p_4,p_5,q_0,q_3,q_4,q_5$. This implies that the (complex-valued) linear form $Q = (q_1{-}iq_2)(\theta^1{+}i\theta^2)$ and the (complex-valued) cubic form $P = (p_1{-}ip_2)(\theta^1{+}i\theta^2)^3$ are well-defined on the $2$-sphere $S^2 = \mathrm{SO}(3)/\mathrm{SO}(2)$, where the $\mathrm{SO}(2)$-subgroup of $\mathrm{SO}(3)$ is the one that is an integral of the forms $\theta^1$ and $\theta^2$. Thus, $P$ and $Q$ are sections of natural complex line bundles over the $2$-sphere.

One then finds that the integral manifolds of the differential system must satisfy $$ \begin{align} v_5 &= -\sinh(\rho)\ (p_4{+}q_5)-\cosh(\rho)\ (p_5{+}q_4)\\ v_6 &= (\cosh^2(\rho){-}\cosh(\rho)\sinh(\rho){-}2)\ p_0 +(\cosh^2(\rho){+}\cosh(\rho)\sinh(\rho){-}1)\ (\cosh(\rho)\sinh(\rho){-}\cosh^2(\rho){+}1)\ p_3\\ &\quad +(\cosh(\rho)\sinh(\rho){-}\cosh^2(\rho){+}{\tfrac{1}{3}})\ q_0 +(\cosh^2(\rho){-}\cosh(\rho)\sinh(\rho))\ q_3\\ v_7 &= (\cosh^2(\rho){-}\cosh(\rho)\sinh(\rho){-}2)\ p_0 +(\cosh^2(\rho){-}\cosh(\rho)\sinh(\rho){+}1)\ (\cosh(\rho)\sinh(\rho){+}\cosh^2(\rho){-}1)\ p_3\\ &\quad +(\cosh(\rho)\sinh(\rho){+}\cosh^2(\rho){-}{\tfrac{1}{3}})\ q_0 +(\cosh^2(\rho){+}\cosh(\rho)\sinh(\rho))\ q_3\ .\\ \end{align} $$

Finally, once these are substituted into the $2$-forms, one finds that the differential system reduces to a pair of second order linear differential equations on $S^2$. Using notation that can be found in a paper of mine (TAMS 290 (1985), 259–271), this second order equation for the tensors $P$ and $Q$ can be written in the form $$ Y(Y(P)) = \tfrac12\ Q + \tfrac23\ X(Y(Q)) - \tfrac13\ X(X({\overline{Q}})), $$ where $X$ and $Y$ are certain first-order differential operators that are invariant under $\mathrm{SO}(3)$ and that generalize $\partial$ and $\bar\partial$ to symmetric $(1,0)$-forms of arbitrary degree. This is an elliptic linear equation of second order that is underdetermined (it is $2$ equations for the $4$ components of $P$ and $Q$). It is invariant under the rotations of $S^2$ and equivalent to the original (overdetermined) EDS on $\mathrm{SO}(3)\times\mathbb{R}$.

It is conceivable that this equation admits an explicit solution in terms of a potential (which would be a section of a rank $2$ vector bundle over $S^2$), but I have not tried to check whether this is true or not.

Added Comment: I have now checked about the possibility of a potential and, miraculously, it turns out that there is a potential: One can show that the solutions of the above equation for $P$ and $Q$ are expressible in the form $$ \begin{align} P &= X(X(L))\\ Q &= L + (XY{+}YX)L + X(X({\overline{L}})\bigr)\\ \end{align} $$ where $L=(L_1{+}iL_2)(\theta^1{+}i\theta^2)$ is an arbitrary complex-valued $(1,0)$-form on $S^2$. Thus, components of $L$ are the two arbitrary functions of $2$ variables predicted by the theory for the general solution.

$L$ is not quite unique. It turns out that one can also add an expression of the form $X(a_1+ib_2)$ to $L$ where $a_1$ is the restriction to the $2$-sphere of a linear function in $\mathbb{R}^3$ and $b_2$ is the restriction to the $2$-sphere of a harmonic homogeneous quadratic function in $\mathbb{R}^3$. This describes the (local and global) ambiguity in the potential $L$ completely.

Thus, with this construction and the above formulae, we have the complete description of the (local and global) integral manifolds of the original EDS in terms of the second and third derivatives of $L$ (which is arbitrary).