Solution to second order differential equation I'm reading a paper in which the authors solve the following equation:
$\frac{d^{2}}{dz^{2}}\hat{p}$($\bf{q}$$,z)$-$q^{2}\hat{p}$($\bf{q}$$,z)$-$\frac{iq_{y}}{(2\pi)^{2}}\delta(z-z_{2})$=0
here $\bf{q}$=$(q_{x},q_{y})$ and $q^{2}$=$q_{x}^{2}$+$q_{y}^{2}$ 
$\hat{p}$($\bf{q}$$,z)$ is the fourier transform of the real function $p$($\bf{s}$$,z)$:
$\hat{p}$($\bf{q}$$,z)$= $\int$$p$$(\bf{s}$$,z)$$e^{-i\bf{q}\cdot\bf{s}}$$d\bf{s}$
$z_{2}$ is a parameter in the differential equatin , $i$ is the imaginary unit.
the autors report the following solution:
$N_{1} cosh(qz)$+$N_{2} sinh(qz)$-$\frac{iq_{y}}{8\pi^{2}q}$$e^{|z-z_{2}|}$
$N_{1}$ and $N_{2}$ are coefficients that depend on the boundary conditions.
when i'm solving that equation with mathematica i'm getting different results:
I'm considering $p$ as being just a function of $z$ and treating $\bf{q}$ as a parameter
 DSolve[p''[z] - (q^2) p[z] - (I*qy/(2 Pi)^2) DiracDelta[z - z2] == 0, 
  p[z], z]

and i get:
p[z] = E^(q z) C[1] + E^(-q z) C[2] + (
 i E^(q z - q z2) qy HeavisideTheta[z - z2])/(8 Pi^2 q) - (
 i E^(-q z + q z2) qy HeavisideTheta[z - z2])/(8 Pi^2 q)

which is different from the one the authors are reporting in the paper.
Please can someone explain me how to corectly solve that equation? 
thanks in advance
 A: First, note that the first two terms of their solution are equivalent to your first two terms under the identification $(C_1,C_2)=\frac{1}{2}(N_1+N_2,N_1-N_2)$. So the only potential problem is in the second halves of the two results; even there, we can recognize a common factor of $-i q_y/8\pi q^2$, so the only issue is the exponentials themselves.
To see the issue more clearly, let's explicitly distinguish the $z<z_2$ and $z>z_2$ cases. We have $$e^{-q|z-z_2|}=\begin{cases} e^{q(z-z_2)}, & z<z_2 \\ e^{q(z_2-z)}, & z>z_2 \\  \end{cases}$$

This can be written in terms of Heaviside step functions as 
$\Theta(z-z_2) e^{q(z_2-z)}+\Theta(z_2-z) e^{q(z-z_2)}$, so that the chice of which sign is in the exponent flips at $z=z_2$. This is nearly the same as the Mathematica answer, except that the step functions in the Mathematica answer are identical. Consequently the Mathematica answer has the second half of the answer being zero for $z < z_2$ and equal to the wrong value for $z > z_2$. So I'd trust the paper's result more than Mathematica...
If you really want Mathematica to do it right, you probably have to introduce the Dirac delta via boundary conditions rather than in the function itself. That is, you should integrate your ODE from $z=z_2-\epsilon$ to $z=z_2+\epsilon$. This will imply a 'jump' condition for the first derivative, which can be used as an internal boundary condition for the problem.
