Explicitly verifying a scattering theory identity I have recently studied scattering theory on a formal level and I think I understand the subject quite well by now. However what I often struggle with is to translate the abstract identities into explicit representations and solve problems with it. I have condensed this issue down to the following example problem, which requires a little bit of algebra but is rather instructive in my opinion. I have provided most of the formulae already, there is probably a conceptual mistake somewhere though.

Consider the one-dimensional Schrödinger equation
$$\left(-\frac{1}{2}\frac{d^2}{dx^2} + V(x)\right)\psi(x) = E\psi(x)$$
with the finite square well potential that is terminated by an infinite barrier on one side
$$ V(x) = \begin{cases}
    \infty, & \text{for } x \leq -L \\
    V_0, & \text{for } -L \leq x \leq 0 \\
    0, & \text{for } 0 \leq x.
  \end{cases}$$
For simplicity assume $V_0<0$.
One set of scattering states for this problem is easily found as
$$\psi^{(+)}(E,x) = \frac{1}{\sqrt{2\pi}}\begin{cases}
    \frac{I(E)\beta}{\alpha}\sin\left(\alpha(1+\frac{x}{L})\right), & \text{for } -L \leq x \leq 0 \\
    e^{-i\sqrt{2E}x} + S(k) e^{i\sqrt{2E}x}, & \text{for } 0 \leq x.
\end{cases}$$
Here, $S(k)$ is the scattering matrix (just a number since there is only the reflection channel here)
$$S(E)=-\frac{\alpha\cot(\alpha) + i\beta}{\alpha \cot(\alpha) - i\beta}$$
and the remaining coefficients are
$$I(E) = - \frac{2i\alpha}{\alpha \cos(\alpha) - i\beta \sin(\alpha)},$$
$$\alpha = \sqrt{\beta^2-2V_0L^2},$$
$$\beta = L \sqrt{2E}.$$
So far so good. Now from formal scattering theory we know that there is also a T-matrix defined by (see e.g. Eq. (7.40) in Newton's book (available on Springer Link))
$$T(E) = \langle\psi_0(E)|V|\psi^{(+)}(E)\rangle.$$
Importantly, the T-matrix is related to the scattering matrix, which in the single channel case takes the simple form (see Eq. (7.58) in Newton's book)
$$S(E) = 1 - 2\pi i T(E).$$
Here, $\psi_0$ is an eigenstate of the free Hamiltonian (i.e. with $V=0$), in our example with the boundary condition at $x=-L$ we get
$$\psi_0(E,x)=\sqrt{\frac{2}{\pi}}\sin\left(\beta(1+\frac{x}{L})\right).$$
Now for our example the overlap integral for the T-matrix can be evaluated in the position representation
$$T(E) = V_0 \int_{-L}^{0} dx \psi_0(E,x) \psi^{(+)}(E,x)$$
and we can plug in our formulae for that. However when substituting the result into the relation to the scattering matrix, it does not hold. I have checked this using Mathematica and manual calculation.

I am clearly doing something wrong. But what? My suspicion is that I have plugged in the wrong states, but I don't know what the right ones are.

EDIT: Following the discussion with TwoBs, here is some more insight on which states should be used. As far as I understand $\psi_0(E,x)$ can just to be an eigenstate of the free Hamiltonian; $\psi^{(+)}(E,x)$ is an eigenstate of the full Hamiltonian but also defined uniquely by the Lippmann-Schwinger equation:
$$|\psi^{(+)}(E)\rangle = |\psi_0(E)\rangle + G^{(+)}(E) V |\psi^{(+)}(E)\rangle,$$
with $G^{(+)}(E) = \frac{1}{E-H_0 + i0^+}$.
The explicit formula I gave for $|\psi^{(+)}(E)\rangle$ above was just some eigenstate of the full Hamiltonian, so the mistake is probably that it does not fulfill the Lippmann-Schwinger equation with the $|\psi_0(E)\rangle$ I used. But which state does? 
 A: This is mostly correct. To make it work, a few inconsistencies in the
definitions and notation need to be taken care of.
First, you are defining the Hamiltonian with $V_0=0$ to be the no scattering
situation. With no scattering, $T=0$, so the $S=1+2i\pi T$ equation was derived
with the convention that $S = 1$ when there is no scattering. Instead
with $V_0=0$, your S contains the phase shift of the wave traveling
from $x=0$ to $x=-L$ and back along with a sign change from the reflection.
To be consistent, you should define $S$ without this extra phase shift:
\begin{equation}
\psi^{(+)}(E,x) = \frac{1}{\sqrt{2\pi}} \left \{
\begin{array}{cc}
\frac{I(E)\beta}{\alpha}\sin\left (\alpha(1+\frac{x}{L})\right),
& -L \leq x \leq 0\\
e^{-i\sqrt{2E}x}-S(E)e^{i\sqrt{2E}2L}e^{i\sqrt{2E}x}
\end{array}
\right .
\end{equation}
with
\begin{equation}
S(E) = \frac{\alpha\cot\alpha+i\beta}{\alpha\cos\alpha-i\beta}
e^{-i\sqrt{2E}2L}
\end{equation}
The Green's function for the $\psi^+(E,x)$ state that satisfies the
Lippmann-Schwinger equation contributes only outgoing waves at infinity.
So this tells you that for large $x$ it should be in the form
$\psi_0(E,x)+Ce^{i\sqrt{2E}x}$, where the $\psi_0(E,x)$ contributes
all of the incoming waves. Your scattering
solution has the incoming wave $e^{-i\sqrt{2E}x}$ multiplied by
$\sqrt\frac{1}{2\pi}$, but in your $\psi_0(E,x)$ this term is
multiplied by $-\frac{i}{\sqrt{2\pi}}e^{-i\sqrt{2E}L}$. So the simplest workaround is
to redefine $\psi_0(E,x)$ to be
\begin{equation}
\psi_0(E,x) = i\sqrt{\frac{2}{\pi}}e^{-i\sqrt{2E}L}\sin\left
(\beta(1+\frac{x}{L})\right) \,,
\end{equation}
so that it is consistent with $\psi^{(+)}(E,x)$ and the Lippmann-Schwinger
equation as TwoBs suggested.
Finally, the on-shell $T$ matrix $T(E)$ should be defined from
the half on-shell $T$ matrix
$T(E',E)=\langle \psi_0(E')|V|\psi^{(+)}(E)\rangle$ by
$\int dk' \delta(\frac{k'^2}{2}-E)T(\frac{k'2}{2},E)$, where
the integral is given by how you form the completeness relation with
your normalization.
Therefore $T(E)= \frac{1}{\sqrt{2E}} \langle \psi_0(E')|V|\psi^{(+)}(E)\rangle$
See, for example S. Weinberg, Quantum Theory of Fields, Cambridge
University Pres, 1995, Eq. 3.2.7.
With these changes, integrating
\begin{equation}T(E)=\frac{V_0}{\sqrt{2E}}\int_{-L}^0 \psi_0(E,x)\psi^{(+)}(E,x)\end{equation}
will give $S(E)=1-2\pi i T(E)$.
Appendix
Here is a more detailed explanation of half on-shell comment above:
If you start with
\begin{equation}
(H_0+V)|\psi^\pm_\alpha\rangle=E_\alpha|\psi^\pm_\alpha\rangle\,,
\end{equation}
where $\alpha$ labels the particular eigenstate $|\phi_\alpha$
and solve
with the boundary conditions that all incoming (+) or outgoing (-)
waves come from the
initial state $|\phi_\alpha\rangle$, you have
\begin{equation}
|\psi^\pm_\alpha\rangle = |\phi_\alpha\rangle
+ \frac{1}{E_\alpha-H_0+i0^\pm} V|\psi^\pm_\alpha\rangle\,.
\end{equation}
Here $|\phi_\alpha\rangle$ is an eigenstate of $H_0$ with energy $E_\alpha$.
If you insert a complete set of such eigenstates of $H_0$, you define the
half on-shell T-matrix
\begin{equation}
T_{\beta\alpha}=
\langle \phi_\beta|\frac{1}{E_\alpha-E_\beta+i0^+} V|\psi^+_\alpha\rangle\,.
\end{equation}
It is half on-shell because $E_\alpha$ is the eigenstate energy, but
$E_\beta$ is not. You can  also solve this where $E_\alpha$ is replaced
with an arbitrary complex number, so that neither energy is
on-shell. The S-matrix is given by
\begin{equation}
S_{\beta\alpha} = \delta(\alpha-\beta)-i2\pi \delta(E_\alpha-E_\beta)
T_{\beta\alpha}\,.
\end{equation}
With the delta function, only terms with
both energies equal contribute,
so these give the on-shell T-matrix. Typically, you solve
for the all the components of the half on-shell T-matrix, but its on-shell components
give you the scattering.
The above delta function is the one I referred to above. It is
the delta function you would have in Fermi's golden rule if
you were to calculate the transition rate (in lowest order Fermi's
golden rule uses the potential matrix elements, but in higher order
these are replaced by the T-matrix elements). The cross section is
the transition rate divided by the incoming flux and you would get
the same sort of relation.
