Ohhh boyy this one was a doozy. 4 years and almost 2 degrees later I can finally answer this question. I'm reviving this ancient thread because it is the first example in google when you search "simple symplectic reduction", and also there are very few examples of applied symplectic reduction out in the wild, so hopefully somebody might find my struggles beneficial. In my original post, I was fine up until Equations 13-15. That is, Equation 12 has been correctly reduced, but there's much more going on under the hood than I initially realized. In a sense, I had successfully applied symplectic reduction, but was not thinking about the resulting BVP correctly. Over the past several years, my understanding of the result has greatly changed.
In my second post on this page, I mentioned I had some difficulties completely decoupling the constants from the original system. That's because, save for some special conditions on the boundary, it is impossible to completely decouple the constants-of-motion from the dynamics. Fundamentally, we're solving a Hamiltonian BVP to indirectly solve the original functional minimization problem. Constants-of-motion have either real physical, or mathematical meaning in context of the original problem and therefore cannot be completely ignored. As a concrete example, say we're in a spacecraft performing an orbit raising manuever. The adjoint of true anomaly is a constant of motion, however the point along the orbit where spacecraft starts and finishes its maneuver is fairly significant in context of the mission and cannot be completely removed.
Are we out of luck? Nope! In trying to discover a compact procedure for application to some real problems and squeezing some benefits out of it, I came up with the idea to separate quantities into two categories: "dynamical" and "non-dynamical" terms. Terms that are dynamical affect the dynamical evolution of a system. Terms that are non-dynamical... well... do not affect the dynamics. Starting off with the original problem:
$$
\begin{aligned}
\min J &= \frac{1}{2} \int_{0}^{t_f} u^2 dt \\
\text{Subject to} : \dot{x}_1 &= x_2 \\
\dot{x}_2 &= u \\
x_1(t_f) &= x_{1f}
\end{aligned}
$$
Plus some other boundary conditions, but I'm only including 1 to illustrate a point and make things too cluttered. Adjoining these to the problem, we get the equivalent functional:
$$
\min J = \frac{1}{2} \int_{0}^{t_f} u^2 + \lambda_1(\dot{x}_1 - x_2) + \lambda_2(\dot{x}_2 - u) dt + \nu(x_1(t_f) - x_{1f})
$$
3 Lagrange multipliers were added: ##\lambda_1##, ##\lambda_2##, and ##\nu##. ##\lambda_1## and ##\lambda_2## are dynamical, whereas ##\nu## is non-dynamical. This is because the ##\lambda##'s will show up in the dynamics but ##\nu## will not. Perturbations in ##\nu## will not affect how the system evolves over time. Following my original procedure in the first post, our dynamical system is:
$$
\begin{aligned}
\dot{x}_1 &= x_2 \\
\dot{x}_2 &= u \\
\dot{\lambda}_1 &= 0 \\
\dot{\lambda}_2 &= -\lambda_1
\end{aligned}
$$
So the ##x##'s and ##\lambda##'s, at this point, are the dynamical states for our system. Hold up, ##x_1## doesn't appear anywhere in the dynamics, so wouldn't that be non-dynamical? And ##\lambda_1## is a constant-of-motion, so is it even a state? It's not a coincidence these facts come in two's, and it's also not a coincidence these states are "adjoint" to one-another, i.e. we don't notice the same thing with ##x_i## and ##\lambda_j## for ##i \neq j##. This all agrees with Noether's theorem. Symplectic reduction requires we get rid of both ##x_1## and ##\lambda_1##. Using ##q_1 = x_1## and ##C_1 = \lambda_1## so we come to my original Eq. 12
$$
H = C_1 x_2 + \lambda_2 u + \frac{1}{2}u^2
$$
Which generates the following system
$$
\begin{aligned}
\text{non-dynamical states} : \; \dot{q}_1 &= x_2 \\ \hline
\text{dynamical states} : \; \dot{x}_2 &= u \\
\dot{\lambda}_2 &= -C_1
\end{aligned}
$$
##C_1## is a dynamical parameter in the sense that small perturbations affect how the system evolves, therefore causing changes in ##x_2## and ##\lambda_2##, which then causes changes in ##q_1##. ##q_1## is non-dynamical because, small perturbations to ##q_1## do not affect anything, other than itself. ##q_1## is a symmetry for the system, and in fact all non-dynamical states are symmetries. Not all symmetries are non-dynamical states, however. In other work, people allude to the fact that ##\dot{q}_1## is solved via quadrature rules, but I've never seen an actual implementation of this. By restricting my definition of a non-dynamical state to be a subset of a symmetries, I was able to write my own general purpose BVP solvers https://github.com/Rapid-Design-of-Systems-Laboratory/beluga/tree/master/beluga/bvpsol for the specific intent of solving Hamiltonian BVPs which have been reduced symplectically. Examples "R2", "R8", and "R18" within the https://github.com/Rapid-Design-of-Systems-Laboratory/beluga/tree/master/beluga/bvpsol/tests are real reduced problems and though they are not symplectic, the resulting BVP structure is identical to symplectic BVPs. To anyone who digs into the code, you'll see that for the BVP solver to take full advantage of our structure here, we need to specifically tell it what is dynamical and what is non-dynamical. I'm publishing a paper on these BVP solvers in ~2 months and will edit this post once it's out.
For those who want to dig deeper into the mathematics behind this, in my framework, trajectories evolve on a principal ##G##-bundle. I'm confident in the structure of my BVP solvers, but rephrasing problems to make use of it remains a challenge. I believe Lie-Poisson mechanics is the correct setting to fully exploit the Hamiltonian system and pose it as a BVP with dynamical and non-dynamical states, but this is ongoing.