Finding approximate eigenfunctions solutions with small eigenvalues This question is about an appendix to chapter 7 of Aspects of Symmetry Erice lectures by Sidney Coleman. We have a TISE for a 1-dimensional simple harmonic oscillator with $\omega = 1$, describing the bottom of a well of a symmetric double potential well, with vanishing minima centered at $\pm a$. At the right minimum, for $a>x\geq 0$,
$$-\frac{1}{2} \hbar \, \partial_x^2 \psi \, +\frac{1}{2}(x-a)^2\psi = E\psi .\tag{A.2.15}$$
For the ground state, we know $E=\frac{\hbar}{2}$ and the even solution $\psi_1 = e^{-(x-a)^2/2\hbar}$. Because of linearity and reflection symmetry, we know there is another solution, $\phi_1$, odd and increasing for which we don't know an analytic form, but for small $\hbar$ its asymptotic approximation is $\phi_1 = e^{(x-a)^2/2\hbar}/(a-x)$. Here he notes that with the normalizations chosen for $\psi_1$ and $\phi_1$, the Wronskian of the two solutions is
$$\phi_1 \partial_x \psi_1 - \psi_1 \partial_x \phi_1 = 2/\hbar.$$
Now, he defines a small dimensionless quantity $\epsilon$ to represent the first order correction to the ground state energy level, by $E=\hbar \, \big(\frac{1}{2} + \epsilon\big)$. Starting from $\psi_1$ and $\phi_1$, solutions to the SE for $\epsilon=0$, we want to find an approximate solution $\psi$ for small $\epsilon$, neglecting $\epsilon^2$ terms.
Here, Coleman says that this can be achieved by a standard method, namely turning the above TISE into an integral equation and iterating once. He immediately writes this solution to be
$$\psi = \psi_1 - \epsilon \int_{x}^{\infty} \mathrm{d}x' \, \psi_1(x')\big[\psi_1(x')\phi_1(x)-\phi_1(x')\psi_1(x)\big].\tag{A.2.26}$$
How did he derive this? What is this "standard method of integrating and iterating once" he's talking about? I substituted this solution $\psi$ back into the SE and after computing the derivatives, using the Wronskian written above and the fact that we know $\psi_1$,$\phi_1$ to be solutions for $\epsilon=0$, I was able to verify by inspection that it is indeed the case that this $\psi$ solves the SE, except for terms of order $\epsilon^2$, which we neglect. But I was not be able to derive this $\psi$  from scratch.
How did he do it?
 A: OK, here is an arguably useful, very schematic trail map of the landscape. I simply lack the patience and gumption to fuss the details... totally irresponsibly cavalier and insouciant of normalizations, signs, etc, with $\hbar=1$, etc, here goes:
$$
(\partial_x^2 -V(x)) \psi(x) = \epsilon \psi (x), 
$$
where I've multiplied by 2, and flipped signs, and shifted $x\mapsto x+a$, etc, and taken the rhs -1 to the left to incorporate into the harmonic potential, etc...
The solutions of the homogeneous equation $\epsilon =0$ are the $\psi_1, ~~\phi_1$, as sketched. The standard conversion of the TISE to a Fredholm integral equation (Mathews & Walker Ch 11) amounts to finding the Green's function to the thus shifted Hamiltonian,
$$
(\partial_x^2 -V(x)) K(x,y)=\delta (x-y),
$$
so  operating  both sides of the above with it and integrating by parts yields the standard Fredholm equation
$$
\int\!\! dy K(x,y)(\partial_y^2 -V(y))(\psi(y)-\psi_1(y)) = \epsilon \int\!\! dy ~K(x,y) \psi(y), \implies \\
\psi(x) = \psi_1(x)+\epsilon \int\!\! dy ~K(x,y) \psi(y)
$$
where the homogeneous solution $\psi_1$ is chosen as a "BC" of the Green's function, as it behaves well (attenuates) for large x Coleman has chosen to consider. (You have seen the analog for the Green's function of the free laplacian in the Born approximation scattering formula.)
The messy technical part is huffing and puffing$^\natural$ to produce something like the Volterra kernel $K(x,y)\sim \theta (y-x) (-e^{(-y^2 +x^2)/2}/x+ e^{(y^2-x^2)/2}/y)$, that  is the Green's function for the oscillator right branch. (It is not the Mehler kernel, nor should it be, given these BCs, but it  vaguely evokes it.) Coleman rarely missed an opportunity to be educational, so if there were obvious shortcuts to this, he would not have resisted sticking a dandy footnote in. But, as you see from Brysk, it is not self-evident.
In any case,  iterating it once (plugging the l.h.s. into the r.h.s.) just yields the trivial term of the Neumann series,
$$
\psi(x) = \psi_1(x)+\epsilon \int\!\! dy ~K(x,y) \psi_1(y) ~~~+O(\epsilon^2).$$

$^\natural$Henry Brysk (1963). "Determinantal Solution of the Fredholm Equation with Green's Function Kernel", J Math Phys 4 (1963) 1536 
