You did not include the proper infinitesimal imaginary part for frequency in the first equation for the Green's function, which would have given the correct pole structures. So it is not clear whether this is time-ordered, retarded or advanced. The retarded Green's function is
$G_R(\omega,p)=\dfrac{\omega+\xi}{(\omega+i\delta)^2-\xi^2-\Delta^2}$
You can either obtain this directly, or perform an analytical continuation from Matsubara(imaginary-time) Green's function.
The imaginary part is
$\Im G=-{\pi(\omega+\xi)}\Big(\delta(\omega-\sqrt{\xi^2+\Delta^2})+\delta(\omega+\sqrt{\xi^2+\Delta^2}) \Big)$
Another problem is that in a superconductor the Green's function is really a matrix. The one you have is probably just $\langle c c^\dagger\rangle$. The density of state is given by the imaginary part of the trace of the matrix Green's function. The full Green's function is given by
$\dfrac{\omega+\xi\tau_z+\Delta\tau_x}{\omega^2-\xi^2-\Delta^2}$
Here $\tau$ are Pauli matrices in the Nambu (particle-hole) space.
For $p$-wave superconductors, well, once you write down the Green's function in momentum space, you have assumed translation invariance and a constant superconducting gap everywhere. Zero-energy state only appears at defects, where the order parameter vanishes (like the core of a vortex, or the edge). So you don't see the zero-energy state from this particular Green's function because there are not any. You have to solve BdG equation in real space, or solve the full Gor'kov equation for the Green's function in real space to obtain the zero energy state.
To take things in order
How is the third Green's function derived? I ask for detailed derivation process.
To reduce clutter and focus on ideas, I will drop a lot of subscripts and such. Note: the clutter I'm removing (e.g. $k,\alpha,\beta$) is still there, I'm just not writing it. I will also drop $\mu$ since I usually don't see it appearing the way you have it (see: https://en.wikipedia.org/wiki/Green%27s_function_(many-body_theory)#Zero-temperature_limit) and it can always be handled by absorbing it into the Fermi function or in shifting the energies. Begin with the second equation
$$G(\omega)=\sum_n \left\{\frac{\langle \psi|c|n\rangle \langle n|c^{\dagger}|\psi\rangle }{\omega-(E_n-E_m)+i0^+}+\frac{\langle \psi|c^\dagger|n\rangle \langle n|c|\psi\rangle }{\omega+(E_n-E_m)-i0^+}\right\}$$
Insert $1=\int d\omega'\delta(\omega'-E_n+E_m)$ to get
$$
G(\omega)=\int d\omega'\sum_n \delta(\omega'-E_n+E_m)\left\{\frac{\langle \psi|c|n\rangle \langle n|c^{\dagger}|\psi\rangle }{\omega-(E_n-E_m)+i0^+}+\frac{\langle \psi|c^\dagger|n\rangle \langle n|c|\psi\rangle }{\omega+(E_n-E_m)-i0^+}\right\}
$$
The delta function enforces $\omega'=E_n-E_m$ so we can substitute this into the expression. Note also that now the numerators in the integral are $A(\omega')$ and $B(\omega')$ respectively. This gives
$$
G(\omega)=\sum_n \left\{\int d\omega'\frac{A(\omega')}{\omega-\omega'+i0^+}+\int d\omega'\frac{B(\omega')}{\omega+\omega'-i0^+}\right\}
$$
Which is precisely your third expression.
Why is the imaginary part of the Green function the spectral function?
Well let's take the imaginary part of $G$.
$$\Im G(\omega) = \Im \int d\omega'\sum_n \left\{\frac{A(\omega')}{\omega-\omega'+i0^+}+\frac{B(\omega')}{\omega+\omega'-i0^+}\right\}
$$
For compactness let's drop the sum on $n$ and the $B(\omega')$ term - they all work the same way. To make the next part clearer, let's also replace $0^+$ by $\lim\limits_{\varepsilon\to0^+}\varepsilon$. So we are trying to analyze
$$
\Im \lim\limits_{\varepsilon\to0^+}\int d\omega'\frac{A(\omega')}{\omega-\omega'+i\varepsilon}
$$
Multiply the numerator and denominator the by $\omega-\omega'-i\varepsilon$ to get
$$
\Im \lim\limits_{\varepsilon\to0^+}\int d\omega'\frac{A(\omega')(\omega-\omega'-i\varepsilon)}{(\omega-\omega')^2+\varepsilon^2}=
$$
$$\Im \lim\limits_{\varepsilon\to0^+}\int d\omega'A(\omega')\left[\frac{-i\varepsilon}{(\omega-\omega')^2+\varepsilon^2}+\frac{\omega-\omega'}{(\omega-\omega')^2+\varepsilon^2}\right]
$$
(Note: taking the limit here gives a version of a Kramers–Kronig relation you ask about)
$A(\omega)$ is real and non-negative for Fermions (exercise left for reader) so the second term has no imaginary part and can be dropped. The first term is $-i\pi$ times a nascent delta function (specifically the Poisson kernel). This allows us to rewrite the above as
$$
\Im \lim\limits_{\varepsilon\to0^+}\int d\omega'A(\omega')\frac{-i\varepsilon}{(\omega-\omega')^2+\varepsilon^2}=\Im -i\pi\int d\omega'A(\omega')\delta(\omega-\omega')=-\pi A(\omega)
$$
Which gives us the known relation $-\frac{1}{\pi}\Im G(\omega,\vec{k})=A(\omega,\vec{k})$ (other conventions may give $-2$ instead of $-1/\pi$). I've brought back the $\vec{k}$ to emphasize that this doesn't hold for $G(\omega)=\sum_{\vec{k}} G(\omega,\vec{k})$ and $A(\omega)=\sum_{\vec{k}} A(\omega,\vec{k})$
The spectral function is for what?
The spectral function gives information about the excitations (particles) in your system. As you can see from the form of the spectral function $A(\omega,k)=\sum_n \delta(\omega-E_n+E_m)\langle \psi|c_k|n\rangle \langle n|c_k^{\dagger}|\psi\rangle $ the spectral function will be peaked at the energies and wavevectors associated with the species in your system. As we showed above, the spectral function allows us to get the Green's function. It can be used to get the filling of the system and information about the density of states. (Note that this applies to noninteracting systems which is what you've presented here. The interacting case is more subtle but the idea holds that the spectral function tells you about the excitations in your system).
What is the relationship between the above thing and the Kramers-Kronig relationship?
There are a few ways to understand the Kramers-Kronig relations. One is that given a function that's analytic in UHP and decays at least as fast as $1/|\omega|$ for large $|\omega|$) then the Kramers-Kronig relations allow you to reconstruct the whole function from knowing just the real part or just the imaginary part. This is essentially how we can use the spectral function (the imaginary part of $G$) to get the full $G$! This is why the third equation makes sense and how we show that the spectral function is (up to a constant) the imaginary part of the Green's function. These relations are a special case of the Sokhotski–Plemelj (https://en.wikipedia.org/wiki/Sokhotski%E2%80%93Plemelj_theorem) and is a useful thing to know for anyone working in condensed matter physics.
Hope that clears things up, cheers!
Best Answer
I cannot comment, that's why I post an answer.
I argee with answers above: you don't need Green's function to solve this equation, you just need to separate variables by entering spherical coordinate system like this $\psi(x,\theta,\phi)=R(r)\psi_{ang}(\theta,\phi)$. After doing this you'll get two separate equations on functions $R(r)$ and $\psi_{ang}(\theta,\phi)$: Bessel equation for $R(r)$, see Bessel functions and equation for angles. Solution of second equation is represented by spherical harmonics.
P.S. It is also easier to search for $R(r) = \frac{\chi(x)}{x}$. Than you'll exactly get Bessel equation.