Prove that this differential equation has reflective symmetric about a line Context
I am studying analytic mechanics using  [1]. In the context of orbital mechanics the authors arrive at the following initial-value problem.
\begin{align}
\frac{d^2 u{(\theta)}}{d\theta^2 } 
 +
u{(\theta)}    = 
-\frac{m}{\ell^2   }\,    \frac{\partial{V( u)}}{\partial{u}}  
\end{align}
with conditions
\begin{align}
\left.V(u)\right|_{u=0} = 0;
\\
\left.u(\theta )\right|_{\theta_r} = u_o;
\\ 
\left.\frac{du(\theta)}{\partial\theta}\right|_{\theta_r} = 0.
\end{align}
The authors claim that, the orbit is invariant under reflection about the line of apsides. From [2], "the apsides refer to the farthest (1) and nearest (2) points reached by an orbiting planetary body (1 and 2) with respect to a primary, or host, body (3). The line of apsides is the line connecting positions 1 and 2."
I do not like the proof given in [1]. I have written my own; and, frankly, I like it less still.
Question
Given the initial-value problem above, prove the solution is symmetric about the line of apsides.
My Attempt
Outline of my proof.
Part 1
I know that the homogeneous solution of the governing differential equation are sinusoidal functions [3]. I use a complete orthogonal set as a basis. I write
$$u = \sum_{n=0}^\infty f_n(u)\,\sin(n\,[\theta - \theta_r]) + g_n(r)\,\cos(n\,[\theta - \theta_r])  .$$
The initial condition will only be satisfied if $f_n(u)=0$ for all $n$. Plugging into the originating differential equation, I have that
\begin{align}
\sum_{n=0}^\infty \left[1 -n^2\right]\,g_n(u)\,\cos(n\,[\theta - \theta_r]) 
     = 
-\frac{m}{\ell^2   }\,    \frac{\partial{V( u)}}{\partial{u}}  
\end{align}
Therefore, the solution must satisfy
\begin{align}
  g_0(u)
+
\sum_{n=2}^\infty \left[1 -n^2\right]\,g_n(u)\,\cos(n\,[\theta - \theta_r]) 
     = 
-\frac{m}{\ell^2   }\,    \frac{\partial{V( u)}}{\partial{u}}  .
\end{align}
From orthogonality,  for all $p\geq 2$,
\begin{align}
  g_p(u)  
     = 
-\frac{m}{\ell^2  \left[1 -p^2\right] }\,    
\frac{
\int_0^{2\,\pi} \cos(p\,[\theta - \theta_r]) \,d\theta
}
{
\int_0^{2\,\pi} \cos^2(p\,[\theta - \theta_r]) \,d\theta
}
\,\frac{\partial{V( u)}}{\partial{u}}  =0
\end{align}
Though unnecessary to determine for this problem, I use the initial conditions to write that
\begin{align}
g_0(r)  
&= 
-\frac{m}{\ell^2   }\,    \frac{\partial{V( u)}}{\partial{u}}   ,~\text{and}
\\
g_1(r) 
&= 
u_o + \frac{m}{\ell^2   }\,    \frac{\partial{V( u)}}{\partial{u}}    
  .
\end{align}
Thus,
\begin{align}
 u(r,\theta) 
&=   
-\frac{m}{\ell^2   }\,    \frac{\partial{V( u)}}{\partial{u}}   + \left[   
u_o + \frac{m}{\ell^2   }\,    \frac{\partial{V( u)}}{\partial{u}}    
 \right]  \,\cos( \theta - \theta_r ) .
\end{align}
With the solution in hand, let's see about reflective symmetry about the line of apsides.   This is done in two parts.
Part 2
First, suppose the ordered triple $(a, \theta_r +b,c)$ is any ordered triple on the graph
of
$$z = u(r,\theta) = g_0(r)   +  g_1(r)\,\cos( \theta - \theta_r )   .$$
Then,
$$c =  g_0(a)   +  g_1(a)\,\cos( b)   .$$
I must show that the ordered triple $(a, \theta_r-b,c)$ also satisfies
$$z  = g_0(r)   +  g_1(r)\,\cos( \theta - \theta_r )   .$$
I write
$$g_0(r)   +  g_1(r)\,\cos( \theta_r-b - \theta_r )  = 
g_0(r)   +  g_1(r)\,\cos( -b )  = 
g_0(r)   +  g_1(r)\,\cos( b ) =  c .$$
Second, suppose the ordered triple $(a, \theta_r +\pi +b,c)$ is any ordered triple on the graph
of
$$z = u(r,\theta) = g_0(r)   +  g_1(r)\,\cos( \theta - \theta_r )   .$$
Then,
$$c =  g_0(a)   +  g_1(a)\,\cos(\pi + b)   .$$
I must show that the ordered triple $(a, \theta_r +\pi-b,c)$ also satisfies
$$z  = g_0(r)   +  g_1(r)\,\cos( \theta - \theta_r )   .$$
I write
$$g_0(r)   +  g_1(r)\,\cos( \theta_r+\pi-b - \theta_r )  = 
g_0(r)   +  g_1(r)\,\cos(\pi -b )  =  g_0(r)   +  g_1(r)\,\cos(\pi +b ) =  c .$$
I conclude that the solution of
\begin{align}
\frac{d^2 u{(\theta)}}{d\theta^2 } 
 +
u{(\theta)}    = 
-\frac{m}{\ell^2   }\,    \frac{\partial{V( u)}}{\partial{u}}  
\end{align}
is symmetric about the line of apsides since:
(1) $(a, \theta_r -b ,c)$ lies on the solution's graph when $(a, \theta_r +b ,c)$ does; and
(2) $(a, \theta_r +\pi -b ,c)$ lies on the solution's graph when $(a, \theta_r + \pi+b ,c)$ does.
Discussion
As I said in the context, I do not, overall, like Goldstein's proof. Why this is so is hard to explain. Though I do like one aspect of his approach that my approach lacks. In my approach, I had to first solve the differential equation, and only afterwards was I able to show the reflective symmetry about the line of apsides. In my opinion my approach keeps me in a weak position. To get where I would like to be, I have to move beyond first solving the differential equations in order to see symmetries. It is important to be able to look at the differential equation and under and analyze the symmetries without solving the differential equation first. This is the approach of Goldstein. However, I simply can not follow along with his argument. The question remains open.
References
[1] Classical Mechanics, Goldstein, Poole, and Safko, 3rd Edition, pg 87.
[2] https://en.wikipedia.org/wiki/Apsis
[3] https://en.wikipedia.org/wiki/Ordinary_differential_equation
 A: If you take $x(\theta) = u(\theta)$ and $y(\theta)=u'(\theta)$, then the given differential equation can be re-written as:
$\textbf{x}'(\theta) = \textbf{F}(\textbf{x}(\theta))$ where $\textbf{F}$ is the appropriate smooth multivariate function. Let $\theta_0$ be the coordinate corresponding to the line. Then using the Picard-Lindelof Existence and Uniqueness theorem:
\begin{align*}
\textbf{x}(\theta) = \textbf{x}(\theta_0) +\int_{\theta_0}^{\theta}\textbf{F}(\textbf{x}(\theta))d\theta = \textbf{x}(\theta_0)+ \lim_{n\rightarrow \infty} \int_{\theta_0}^{\theta}\textbf{F}(\textbf{x}_n(\theta))d\theta
\end{align*}
where:
\begin{align*}
\textbf{x}_0(\theta) &= \textbf{x}(\theta_0)\\
\textbf{x}_{n+1}(\theta) &= \textbf{x}(\theta_0) + \int_{\theta_0}^{\theta}\textbf{F}(\textbf{x}_n(\theta))d\theta
\end{align*}
for some neighbourhood of $\theta_0$. All this is to say that for a given $\textbf{x}(\theta_0)$ and $\textbf{F}$, the solution is uniquely determined(for some neighbourhood)! Suppose you get a solution for $[\theta_0-\delta,\theta_0+\delta]$. Consider a new function $\textbf{y}$ as $\textbf{y}(\theta) = \textbf{x}(2\theta_0-\theta)$. Then you can check that $\textbf{y}(\theta_0) = \textbf{x}(\theta_0)$ and that it satisfies the same differential equation for $\theta \in [\theta_0-\delta,\theta_0+\delta]$. By the uniqueness theorem, $\textbf{y}(\theta) = \textbf{x}(\theta)$, and I think this is what the book meant by the assertion.
Edit: All of this only guarantees local existence and uniqueness, and a different variant of the theorem would be required to prove the symmetry holds for the full domain.
