Seeking the most general solution of a first order nonlinear PDE in two independent variables I am seeking the most general form of the solution of the following PDE: $$\left(\frac {\partial G(u,v)}{\partial u}\right)^2+\left(\frac {\partial G(u,v)}{\partial v}\right)^2= F(u)$$where $F(u)$ is  a function of $u$ only. A sufficient solution is $$G(u,v)= e ^{(\pm \alpha u+\beta)}\sin \alpha (v+\gamma)+\epsilon$$ where $\alpha,\beta,\gamma,$ and $\epsilon$ are real constants. 
Furthermore, $F$, $G$, and $u$, $v$ are real-valued, with $F\ne constant$, $0\lt F \lt \infty$, $0\le u \lt \infty $, and $0\le v \le \pi$, where $u$ and $v$ are the oblate spheroidal coordinates like those found in Field Theory Handbook, Second Edition by Moon and Spencer. Is there a more general solution?
The above paragraph from an edit renders irrelevant a lot of the answers below even though the answers are excellent. My apologies for not stating all of the necessary information in the original question.
 A: Robert Bryant alluded to this in a comment but it's maybe worth expanding on.
Your PDE is a type of eikonal equation, which shows up in geometric optics.
You didn't specify the domain or boundary conditions so I'll assume some nice compact domain $\Omega$ with $C^1$ boundary and that $G|_{\partial\Omega} = 0$.
For a particular case I'll take $F = 1$, so the PDE is
$$|\nabla G|^2 = 1$$
although you can generalize this quite a bit.
A weak solution to this problem is
$$G(\mathbf x) = \mathrm{dist}(\mathbf x, \partial\Omega) \equiv \mathrm{min}_{\mathbf y \in \partial\Omega}|\mathbf x - \mathbf y|.$$
Things get more interesting if $F$ (which we can identify with the squared index of refraction of the medium) depends on space, but in any case my point is that the signed distance function to any nice domain can be a solution of your problem.
A: This idea is a little weird and doesn't totally solve your problem but lets see how it goes. We begin by factoring the expression as:
$$ \left( \frac{\partial G}{\partial u} + i \frac{\partial G}{\partial v} \right)   \left( \frac{\partial G}{\partial u} - i \frac{\partial G}{\partial v} \right) = F$$
So this becomes a system of equations
$$  \left( \frac{\partial G}{\partial u} + i \frac{\partial G}{\partial v} \right)   = u_1 \\ \left( \frac{\partial G}{\partial u} - i \frac{\partial G}{\partial v} \right) = u_2 \\ u_1 u_2 =F $$
If my recollection is correct we can solve each individual equation there by considering just ONE solution and then solving the kernel problem. I.E. the general solution to:
$$  \left( \frac{\partial G}{\partial u} + i \frac{\partial G}{\partial v} \right)   = u_1 $$
Is given by finding a test solution $q_1$ and then considering all the solutions of
$$  \left( \frac{\partial G}{\partial u} + i \frac{\partial G}{\partial v} \right)   = 0 $$
Which takes the form $G = H(u+iv)+q_1$ where $H$ is an arbitrary smooth function of your choice. (Or maybe not even smooth?)
For the other equation it would be $G = S(u-iv) + q_2$ where again $S$ is arbitrary.
So now to accept a $H,S$ pair as a solution it must be that our two $G$ expressions are equal in other words:
$$ q_1 - q_2 = S(u-iv) - H(u+iv) $$
This is where things get a little wonky. Depending on your choice of $H,q$ that might not be solvable.
A: This is really not a new answer as it stems out from the nice idea of Sidharth Ghoshal. By looking carefully to his Ansatz, by putting $z=u+iv$ and $\bar{z}=u-iv$, the linear PDEs he obtains can be written as
$$
\begin{cases}
\dfrac{\partial G}{\partial \bar{z}} = u_1 \\
\\ 
\dfrac{\partial G}{\partial {z}} = u_2  
\end{cases}
$$
Considering the case $u_1(u,v)=u_2(u,v)=\delta(u,v)$, by the solution of the $\bar\partial$-problem in one complex variable we have the particular solution
$$
\hat{G}(z,\bar{z})=\frac{1}{\pi z}+\frac{1}{\pi\bar{z}}\label{1}\tag{1}
$$
Thus, a general solution of the original nonlinear first order PDE can be constructed by using \eqref{1} and any couple of factors $(u_1, u_2)$ such that $u_1u_2 = F$ as a sum of convolutions
$$
\begin{split}
G(u,v) & =G(z,\bar{z})=\frac{1}{\pi}\left[\ \int\limits_{\Bbb C} \frac{u_1(\zeta,\bar{\zeta})}{z-\zeta}{\mathrm{d}\bar\zeta\wedge\mathrm{d}{\zeta}} +\int\limits_{\Bbb C} \frac{u_1(\zeta,\bar{\zeta})}{\bar{z}-\bar\zeta}{\mathrm{d}\bar\zeta\wedge\mathrm{d}{\zeta}}\right]\\
&= \frac{1}{\pi}\left[\ \int\limits_{\Bbb R^2} \frac{u_1(u, v)}{(u-x)+i(v-y)}{{\mathrm{d}u\mathrm{d}v}} +\int\limits_{\Bbb R^2} \frac{u_1(u, v)}{(u-x)-i(v-y)}{\mathrm{d}u\mathrm{d}v}\right].
\end{split}\label{2}\tag{2}
$$
where $\zeta=x+iy$, $\bar\zeta=x-iy$.
Notes

*

*In the deduction of formula \eqref{2} we have not used the fact that $F=F(u)$ only as stated in the OP: therefore the formula covers more cases than the one required in the OP.

*Again formula \eqref{2} is written as if $u_1$ and $u_2$ were $L^1(\Bbb R^2)$: however the formula has a meaning also if $u_1$ and $u_2$ are (in some sense) generalized functions, thus the solution is valid even if $u_1$ and $u_2$ are various kind of generalized functions, provided their product is defined and is equal to the given (generalized) function $F$. For example, \eqref{2} has a a well defined meaning even for $u_1, u_2, F \in L^1_\mathrm{loc}(\Bbb R^2)$.

