I have been trying to solve this (see image) equation 1 using the finite difference method with the initial conditions as given in equations (2,3 & 4). Note that the initial condition is z = -2400 nm for all values of the r range. Is it possible to solve this equation (1) with this information? I have been trying without success. Especially, how to approximate r-dependent derivatives without knowing anything about initial conditions, i.e., E_p, E_s at r = -1800 nm for all values of z? Any suggestions, please?

No. You can’t solve a second-order equation in z given only the value of the solution at one z boundary. You need some additional information.

Most likely, since this is a wave equation (a 2d scalar Helmholtz equation), you want to assume outgoing/radiation boundary conditions as |z|,r \to \infty (the r=0 boundary is determined by the cylindrical symmetry). Then you can use a variety of techniques, e.g. finite differences with PML boundaries (as in this example Julia code ), finite elements with PML (as in this Gridap tutorial), or perhaps integral-equation methods (where the outgoing boundaries are implicitly specified via the Green’s function).

@stevengj Thank you for the quick response. I will write more after reading about PML boundary conditions as I am not aware of the theory behind this method. However, just to make a brief response, this equation comes in CARS microscopy [J. Opt. Soc. Am. B 17, 1678 (2000)] where the authors use the integral transform method to the equation (1) i.e., applying Hankel transform to eliminate the r-dependent derivatives in cylindrical coordinates and solving the resulting second-order differential equation in one variable (z) by numerical integration. However, I am trying to solve the same equation (1) by using the finite difference method. Here, the only known initial fields (pump, Stokes) at the front face of the sample cell (assuming a square cell in r, z coordinates) are given by equations 2,3.

I will write more after going through PML theory later. thank you again.

I haven’t looked at this paper, but it sounds like they are effectively using a kind of integral-equation method, which implicitly incorporates outgoing boundary conditions via the Green’s function.

That is right. they use Hankel transform to simplify the radial derivatives by using the known relation of Hankel transform property, that is HT [ d^2/dr^2 + 1/r * d/dr ] = - 4pi rho^2. This way they convert the second-order differential equation in two variables (r , z) into a second-order diff. eq. in single variable, z while rho (conjugate of r) becomes a parameter.