I am attempting to follow this tutorial in the documentation on using FEM to solve PDEs. I am having difficulty understanding how to formulate the Neumann boundary condition for my free-boundary problem. In the notation of the tutorial, it seems integral to be able to write down the equation in the form (I drop the "time" terms since they are not relevant to my question) \begin{equation} \nabla \cdot (-c(t,X,u,\nabla_X u) \nabla u-\alpha (t,X,u,\nabla _Xu) u+\gamma (t,X,u,\nabla _Xu)) + \beta (t,X,u,\nabla _Xu)\cdot \nabla u+a(t,X,u,\nabla _Xu) u - f(t,X,u,\nabla _Xu) = 0\ , \end{equation}

since the term in the first set of parentheses is used to define the Neumann boundary condition:

$$ \overset{\rightharpoonup }{n}\cdot (c \nabla u-\gamma +\alpha u)=g-q u\ . $$

Well, my PDE for $F = F(p, w)$ is

$$ \frac{1}{2}\sigma^2 p^2 \partial_p^2 F + \mu p \partial p + (\rho w + \xi) \partial_w F + (p \partial_w F)^\gamma - \beta F = 0\ , $$

where all the lowercase Greek letters are constants.

(i) I am first of all unclear as to how best to re-write this in the standard form above, since I have a first-order derivative in $w$ but a second-order derivative in $p$. In fact, I can see multiple ways to re-write the equation, but am unsure which one is correct.

(ii) Once this has been accomplished, I would like the following "smooth-fit" conditions to hold on the boundary $\partial \Omega$ of the domain over which I solve the PDE:

$$ F = G,\ \quad \partial_p F = \partial_p G, \quad \partial_w F = \partial_w G. $$

Here $G(p, w)$ is a known function. The first condition is a Dirichlet condition which is straightforward to implement. How would I formulate the latter two conditions in the Mathematica syntax for the Neumann boundary condition?