I'm working on an assignment for my quantum mechanics class and I've arrived at a nonlinear inhomogeneous partial differential equation for a complex function $S:\mathbb{R}^2\to\mathbb{C}~;~S:(x,t)\mapsto S(x,t)$ (This is actually the Quantum Hamilton-Jacobi equation in one spatial dimension) $$\partial _{t} S+\frac{1}{2m} |\partial _{x} S|^{2}+U =\frac{\mathrm{i} \hbar }{2m}\partial ^{2}_{x}S\tag{1}$$ Here, $\mathrm{i}$ is the imaginary unit, $\hbar$ and $m$ are positive real constants (reduced Planck's constant and the mass of the particle), the $|~|$ denotes complex norm, and $U:\mathbb{R}^2\to\mathbb{R}$ is a real function (potential energy). My professor has given me the following instruction:
Expand $S(x,t)$ in terms of the small parameter $\hbar$ $$S(x,t)=S_0(x,t)+\hbar S_1(x,t)+\hbar^2S_2(x,t)+\dots$$ Show that $$\partial_tS_0+\frac{|\partial_x S_0|^2}{2m}+U=0\tag{2}$$ $$\partial_t S_1+\frac{1}{m}(\partial_x S_0)(\partial_x S_1)=\frac{\mathrm{i}\hbar}{2m}\partial_x^2S_1\tag{3}$$
I'm not really sure how he wants us to approach the "expand in terms of the small parameter" part. It looks kind of like a Taylor series, but I don't know how they are related. I've done small parameter expansions on single variable Fredholm integral eigenvalue problems, but I really don't know how I can apply that here. I think my question is kind of similar to this one and this one, but as for the first I'm not sure how to state my equation as an eigenvalue problem and as for the second the only answer on that post is very unclear, and I can't see a way to use geometric series as Ted Shifrin mentioned in the comments. This is from a quantum mechanics class, so we haven't been given a general treatment of the mathematical methods needed for this. "Perturbation theory" is covered later in the course, but only really on linear eigenvalue problems. I have another course that deals with this sort of thing, but its focused mostly on integral equations and on linear single variable differential equations, so I don't know how to apply those methods here.
How should I approach this??
CORRECTIONS:
My professor has issued a revised assignment (still containing errors I might add) with some corrections to equations (1),(2), and (3). They now agree with @Sal 's work. They read $$\partial _{t} S+\frac{1}{2m} (\partial _{x} S)^{2}+U =\frac{\mathrm{i} \hbar }{2m}\partial ^{2}_{x}S\tag{1*}$$ $$\partial_tS_0+\frac{(\partial_x S_0)^2}{2m}+U=0\tag{2*}$$ $$\partial_t S_1+\frac{1}{m}(\partial_x S_0)(\partial_x S_1)=\frac{\mathrm{i}}{2m}\partial_x^2S_1\tag{3*}$$
Obligatory first remark (apologies in advance) $\hbar$ is a dimensionful constant, and so is not 'small' or 'large'. I may choose to work in units where $\hbar=10$, where it is certainly not true that $\hbar^n \to 0$ for large $n$. Let me rewrite the equation, setting $m=1/2$ and $\hbar=1$ for my convenience$^\dagger$
$$ \partial_t S +(\partial_x S)^2 +U =i\partial_{xx}S $$
We wish to solve this equation with the perturbative ansatz
$$ S=\sum_{n=0}^\infty S_n \varepsilon^n $$
That is, we claim the solution may be written as a power series in the small parameter $\varepsilon$: a dimensionless parameter that we have introduced. We wish to solve for the $S_n$ recursively. The next step is to write an $\varepsilon$ in the original equation. Following the spirit of the question, we choose
$$ \partial_t S +(\partial_x S)^2 +U =i\varepsilon \partial_{xx}S $$
The strategy is now straightforward. We substitute our ansatz for $S$ into this equation, and equate powers of $\varepsilon$
$$ \sum_{n=0}^\infty \partial_tS_n \varepsilon^n + \left( \sum_{n=0}^\infty \partial_xS _n \varepsilon^n \right)^2+U=i\sum_{n=0}^\infty \partial_{xx}S_n \varepsilon^{n+1} $$
Products such at the second term on the left are expanded using Cauchy's product, however, we will not need that to get the two terms required in your question. Note the lowest power of $\varepsilon$ on the right is one. The coefficients of $\varepsilon^0$ in the equation are
$$ \partial_t S_0+(\partial_xS_0)^2+U=0 $$
This is your first equation. In principle, you solve this for $S_0$. So presume $S_0$ is known. Look at coefficients of $\varepsilon^1$
$$ \partial_t S_1 + \partial_x S_0 \partial_x S_1 + \partial_x S_1 \partial_x S_0 + 0 = i \partial_{xx}S_0 $$ $$ \partial_t S_1 + 2 ( \partial_x S_0 )(\partial_x S_1) = i \partial_{xx}S_0 $$
Which is your second equation. The middle terms arise from picking the coefficients of $\varepsilon^1$ from the product
$$ (\partial_x S)^2=\left(\sum_{n=0}^\infty \partial_xS _n \varepsilon^n \right)^2 =(\partial_x S_0 + \varepsilon \partial_x S_1 + \varepsilon^2 \partial_x S_2+\dots) (\partial_x S_0 + \varepsilon \partial_x S_1 + \varepsilon^2 \partial_x S_2 +\dots) $$
In principle, since $S_0$ is known, you can now solve the equation for $S_1$, and so on.
$\dagger$ The ansatz $\psi=e^S$ occurs frequently in perturbative analysis. With this ansatz, we get the equation in your post from the Schrodinger equation, without the complex norm. I assume this was a typo, or that what was meant is: $\nabla S \cdot \nabla S$.