Expansion of a plane wave,  Exp[ikz],  
        into spherical Bessel Functions

Regular spherical Bessel functions

j0[x_] := Sin[x]/x
j1[x_] := Sin[x]/x^2 - Cos[x]/x
j2[x_] := (3/x^2 - 1)*Sin[x]/x - 
        3  Cos[x]/x^2
j3[x_] := (15/x^4 - 6/x^2)*Sin[x] + (1/x - 15/x^3)*
          Cos[x]

Singular spherical Bessel functions

\(\(\(
\)\(n0[x_] := Cos[x]/x

    n1[x_] := Cos[x]/x^2 + Sin[x]/x

    n2[x_] := \((3/x^2 - 1)\)*Cos[x]/x + 3  Sin[x]/x^2

    n3[x_] := \(\((15 - 6\ x\^2)\)\ Cos[x] - x\ \((\(-15\) + x\^2)\)\ Sin[x]\)\/x\^4\)\)\)

' Outgoing' spherical Bessel functions

h0[x_] := 
      Exp[I*x]/(I*x)

    h1[x_] := 
      Exp[I*x]*(1/(I*x^2) - 
            1/x)

    h2[x_] := 
      Exp[I*x]*(I/x - 3/(x^2) - 
            3*I/x^3)

    h3[x_] := ((Exp[I*x]* (-15 I - 15 x + 6 I (x^2) + (x^3)))/(x^4))

Spherical harmonics Y^L_  0 (cos (θ))

y0[x_] := 1/Sqrt[4*π]
y1[x_] := Sqrt[3/(4*π)]*x
y2[x_] := 
      Sqrt[5/(16*π)]*(3*x^2 - 1)
y3[x_] := Sqrt[7/π]/4*(5*x^3 - 3*x)

r = 
      Sqrt[z^2 + x^2]

    k = 0.1

    c0 = Sqrt[4*π]

    c1 = 2*I*Sqrt[3*π]

    c2 = -Sqrt[20*π]

    c3 = -2*I*Sqrt[7*π]

    a0 = -Conjugate[c0]*j0'[k*r]/h0'[k*r]

    a1 = -Conjugate[c1]*j1'[k*r]/h1'[k*r]

    a2 = -Conjugate[c2]*j2'[k*r]/h2'[k*r]

    a3 = -Conjugate[c3]*j3'[k*r]/h3'[
            k*r]

r := Sqrt[x^2 + z^2]

Expansion of a plane wave Exp[i*k*z]

f[x_, 
        z_] := (c0*j0[k*Sqrt[x^2 + z^2]])*
          y0[z/Sqrt[x^2 + z^2]] + (c1*j1[k*Sqrt[x^2 + z^2]])*
          y1[z/Sqrt[x^2 + z^2]] + 
+(c2*
                j2[k*Sqrt[x^2 + z^2]])*
          y2[z/Sqrt[x^2 + z^2]] + (c3*j3[k*Sqrt[x^2 + z^2]])*
          y3[z/Sqrt[x^2 + z^2]]

Comparison of absolute values with the exact solution:

Plot[{Abs[f[0, z]], Abs[Exp[I*k*z]]}, {z, -100, 
        100}, PlotPoints -> 30, PlotRange -> {0, 2.3}, 
      AxesLabel -> {z, (ψ^2)}]

(solution OK for  | 
        z | <20)

Plot[{Abs[f[x, 0]], Abs[Exp[I*k*x]]}, {x, -100, 
        100}, PlotPoints -> 30, PlotRange -> {0, 1.1}, 
      AxesLabel -> {x, (ψ^2)}]

(solution OK for  | 
        x | <15)

Solution ψ fulfilling Neumann conditions on the sphere r = 1 for the first four spherical harmonics
(y - direction suppressed)

χ,  for comparison,  
          has an exact ingoing wave and first four harmonics of the outgoing wave

ψ[x_, 
        z_] := (Conjugate[c0]*j0[k*r] + a0*h0[k*r])*
          y0[z/r] + (Conjugate[c1]*j1[k*r] + a1*h1[k*r])*
          y1[z/r] + (Conjugate[c2]*j2[k*r] + a2*h2[k*r])*
          y2[z/r] + (Conjugate[c3]*j3[k*r] + a3*h3[k*r])*
          y3[z/r]
χ[x_, z_] := 
      Exp[-I*k*z] + (a0*h0[k*r])*
          y0[z/r] + (a1*h1[k*r])*y1[z/r] + (a2*h2[k*r])*
          y2[z/r] + (a3*h3[k*r])*y3[z/r]

Phases of ingoing waves:

ContourPlot[
      Arg[Exp[I*k*z]], {x, -50, 50}, {z, -50, 50}, 
      Contours -> 8, PlotPoints -> 20]

ContourPlot[Arg[f[x, z]], {x, -50, 50}, {z, -50, 50}, 
      Contours -> 8, PlotPoints -> 20]

Phases of full solutions:

ContourPlot[Arg[ψ[x, z]], {x, -50, 50}, {z, -50, 50}, 
      Contours -> 8, PlotPoints -> 30, 
      FrameLabel -> {"x", "z"}]

ContourPlot[
      Arg[Conjugate[f[x, z]]], {x, -50, 50}, {z, -50, 50}, 
      Contours -> 8, PlotPoints -> 30, 
      FrameLabel -> {"x", "z"}]

Amplitude of the ingoing wave:

Plot3D[Abs[f[x, z]]^2, {x, -30, 30}, {z, -30, 30}, 
      PlotPoints -> 40, PlotRange -> {0, 3}]

Amplitude of the solution:

Plot3D[Abs[ψ[x, z]]^2, {x, -30, 30}, {z, -30, 30}, 
      PlotPoints -> 40, PlotRange -> {0, 3}]

-x^2*D[n2[x]/x^2, x]

Simplify[%]