Webster equation is a popular generalization of 1D wave equation used for ducts of variable cross-section $S \equiv S(x)$. Assuming harmonicity in time, the spatial equation for propagation of variable $\Phi(x)$ is:
$$ \frac{\partial^2 \Phi}{\partial x^2} + \frac{\partial (\ln S(x))}{\partial x}\frac{\partial \Phi}{\partial x} + k^2\Phi^2 = 0 $$
What FD scheme would be appropriate assuming that we have boundary conditions at $\Phi(0)$ and $\Phi(L)$, i.e. no $\partial_x\Phi(0)$ etc.?
This is not something I've worked on, so I do not know if there is a recommended approach to this, only how I would approach this problem.
I'm assuming you can evaluate $C(x) = (\ln(S(x)))_x$ wherever you want.
Since $\phi$ is prescribed at $0$ and $L$, I would discretise on $x_i = L\frac{i}{N+1}$ for $i$ from $1$ to $N$. Then the finite difference operators could be combined into a matrix $A$ and boundary condition correction term vector $\mathbf{b}$ with $$\begin{gathered} A = \frac{1}{\Delta x^2}\delta^2_x + \frac1{2\Delta x}\mathrm{diag}(C(x_i))\delta_{2x}\,,\\ \mathbf{b}_i = \begin{cases} \left(\frac1{\Delta x^2} - \frac1{2\Delta x}c(x_1)\right)\phi(0)\quad&\text{if}\quad i=1\,,\\ \left(\frac1{\Delta x^2} + \frac1{2\Delta x}c(x_N)\right)\phi(L) &\text{if}\quad i=N\,,\\ 0&\text{otherwise}\,, \end{cases} \end{gathered} $$ so that $$\left[\phi_{xx} + C(x)\phi_x\right]_{x_i} \approx (A\underline{\phi} + \mathbf{b})_i\,.$$
To find the solution, I would use Newton's method, looking for the root of $$F(\underline{\phi}) = A\underline{\phi} + \mathbf{b} + k^2 \mathrm{diag}(\phi_i^2)$$ which has tridiagonal Jacobian $$\mathrm{J}(\underline{\phi}) = A + 2k^2 \mathrm{diag}(\phi_i)\,.$$
(For completion) If there were a Neumann boundary condition, for example $$\phi_x(0) = \alpha\,,$$ I would include $\phi(0)$ in the solution vector $\underline{\phi}$, then I would construct $\mathrm{b}$ by considering a ghost point at $x = -\Delta x$ satisfying $$\frac{\phi(\Delta x) - \phi(-\Delta x)}{2\Delta x} = \alpha\,.$$ I believe mixed boundary conditions wouldn't be much more complicated than non-zero Neumann.
