Choosing initial condition for Hamilton-Jacobi PDE from initial $x$ and $p$ For separable solutions to Hamilton-Jacobi PDE (say in 2D), we treat the Hamilton's principal function $S$ as
$$S= W(x) + W(y) - E*t$$
and treat the separate parts as constants and find $W(x)$, $W(y)$, and finally find $S$, then find the positions or whatever needed by differentiating $W$ with respect to the constant.
But, when the equation is not separable we need to solve the equation numerically and for that initial conditions shall be needed. What are the initial conditions for Hamilton-Jacobi equation, say, for a simple kinetic energy plus potential energy type hamiltonian, and how to find it?
And from the boundary condition, after finding $S$, how can we find the characteristics of the motion, say position as a function of time?
For example, I wish to solve the HJ equation in harmonic oscillator numerically to get position and momentum as functions of time under the initial condition: 
At $$t=0,$$
$$x=1$$
$$p=0$$
Now how shall I proceed numerically, not asking the numerical method of solving the pde, but how to impose the initial conditions on action, and after getting action as a function of position and momentum, how to get $$x$$ and $$p$$ as functions of time.
 A: From Goldstein's book:

[The] technique is to seek a canonical transformation from the coordinates and momenta at time $t$ to a new set of constant quantities, which may be the $2n$ initial values $(q_0,p_0)$ at $t=0$. With such a transformation the equations relating the old and new variables are the desired solution of the mechanical problem: $$q = q(q_0,p_0,t)$$ $$p = p(q_0,p_0,t)$$

Goldstein lays it out very nicely. I will follow his process step-by-step for the harmonic oscillator. We start with the Hamilton-Jacobi equation
$$\frac{1}{2} S^{(1,0)}(x,t)^2+S^{(0,1)}(x,t)+\frac{x(t)^2}{2}=0$$
This has the solution (I used Mathematica)
$$S = \frac{1}{2} \left(x \left(-\sqrt{-x^2+c_1}\right)-c_1 \tan
   ^{-1}\left(\frac{x}{\sqrt{-x^2+c_1}}\right)\right)-\frac{c_1 t}{2}$$
where I have taken the second additive constant of integration to be 0 without loss since only derivatives of $S$ appear in the hamilton-jacobi equation.
The first transformation equation is
$$p = \frac{\partial S}{\partial x}$$
whence
$$p = -\sqrt{-x^2+c_1}$$
the second transformation equation is (I use $\beta$ to emphasize that it is a constant, which is the whole point of this procedure)
$$X = \beta =  \frac{\partial S}{\partial c_1}$$
$$\beta = \frac{1}{2} \left(-t-\tan ^{-1}\left(\frac{x}{\sqrt{-x^2+c_1}}\right)\right).$$
Solving for $x$, we get (there were actually two branches in this whole calculation but I was lazy and only took one, so if you're following very closely that is how the sign ambiguity would be resolved):
$$x = -\sqrt{c_1} \sin (t+2\beta)$$
plug this back into the equation for $p$ and find that
$$p = -\sqrt{c_1}\cos{(t+2\beta)}$$
to fix the constants merely plug in $t = 0$ as well as your initial conditions. The result is
$$c_1\to 1, \qquad \beta \to \frac{1}{2} \left(\frac{\pi }{2}+2 \pi  n\right)$$
A: The Hamilton–Jacobi (HJ) equation is a non-linear 1st order PDE. The flow is typically found using the method of characteristics starting from some initial conditions on a suitable (codimension-1) Cauchy hypersurface. Keep in mind that the flow parameter might not be physical time, and the word "initial" therefore should be interpreted accordingly.
