Assuming typing error in region II (should be -(d/2 +a) < x < -d/2) the potential looks like this:
If it's not what your potential well looks like maybe you would like to post a picture.
The solutions to this problem are obtained by solving one dimensional time independent Schrodinger equation (TISE):
$\hat{H}\Psi(x)=E\Psi(x)$
Where $\hat{H}$ is Hamiltonian, $\Psi$ is an eigen function and $E$ an eigen value. We can further expand the hamiltonian operator
$\hat{H}=\hat{T}(x)+V(x)=\frac{-\hbar^2}{2m}\frac{d^2}{dx^2}+V(x)\,,$
where $\hbar=h/(2\pi)$ is reduced Planck constant.
If you are not sure how to solve the TISE please see finite potential well. If you never did that I strongly recomend that you find it in your text book and first understand the solution of a single finite well. If you don't know what is Schrodinger equation or how the hamiltonian is defined in QM I recomend starting with that. I will describe a general process of solution for the double well. The full solution is rather lengthy because there are many cases. Please let me know which section you don't understand and I'll expand it.
The energy of the particle is $E$. Assuming $V_1<V_2$. We divide the solution into 3 cases:
a)$E<V_1$ (discrete)
b)$V_1<E<V_2$ (discrete)
c)$V_2<E$ (continuous)
For transition probabilities you probably don't need to bother with the discrete spectrum.
In the later I'll use:
$$
V(x)=
\begin{cases}
V_1&x < -(\frac{d}{2} + a)\\
0&-(\frac{d}{2} + a) < x < -\frac{d}{2}\\
V_1& -\frac{d}{2} < x < 0\\
V_2& 0 < x < \frac{d}{2}\\
0&\frac{d}{2} < x < \frac{d}{2} + a\\
V_2&x > \frac{d}{2} + a\\
\end{cases}
$$
We can rewrite the TISE:
$$\frac{-\hbar^2}{2m}\frac{d^2\Psi(x)}{d^2x}+V(x)\Psi(x)=E\Psi(x)$$
$$\frac{d^2\Psi(x)}{d^2x}+\frac{2m(E-V(x))}{\hbar^2}\Psi(x)=0$$
We need to do the folowing part separately for a)$E<V_1$ and b)$V1<E<V2$
Let's take case a)
Let's take a look at region I.
In this region $E<V(x)=V_1$. We define wave vector
$$k_I=\frac{\sqrt{2m(V_1-E)}}{\hbar}$$
(We want $k_I$ to be real number)
We solve the equation:
$$\frac{d^2\Psi_I(x)}{d^2x}-k_I^2\Psi_I=0$$
which leads to solution:
$$\Psi_I=A_Ie^{k_Ix}+B_Ie^{-k_Ix}$$
This applies also for regions IIIa,IIIb,V because in all these regions E
Let's take a look at region II.
In this region $E>V(x)=0$
We define wave vector
$$k_{II}=\frac{\sqrt{2m(E-0)}}{\hbar}$$
We solve the equation:
$$\frac{d^2\Psi_{II}(x)}{d^2x}-k_{II}^2\Psi_{II}=0$$
which leads to solution:
$$\Psi_{II}=A_{II}\sin(x)+B_{II}\cos(x)\,,$$
Similarly for region IV.
Case b) is analogical to a). If $V_1=V_2$ then you only need to calculate case a)
Coefficients
You can calculate the numbers $A_I,B_I,A_{II},B_{II}$ from boundary conditions:
$\Psi_I(-\infty)=0$=>\Psi_I(x)=
$\Psi_I(-(d/2 + a))=\Psi_II(-(d/2 + a))$
$\frac{d\Psi_I}{dx}(-(d/2 + a))=\frac{d\Psi_II}{dx}(-(d/2 + a))$
$\vdots$
$\Psi_V(\infty)=0$
This system of equations has non-zero solution only for certain values of $E$. You can calculate the deterinant in similar manner as for the finite potential well and find values of $E$ for which the determinant is zero. This CAN'T BE SOLVED ANALYTICALLY = you need computer for this step.
Hamiltonian matrix (Overlap)
You can't define ${H=T+V_1+V_2}$ there is no such thing as total hamiltonian. The hamiltonian is:
$$\hat{H}=\frac{-\hbar^2}{2m}\frac{d^2}{dx^2}+V(x)$$
You just plug this in to your integral and get for example:
$$H_{12}=\int_{−(d/2+a)}^{−∞}Ψ_1(x)\frac{-\hbar^2}{2m}\frac{d^2}{dx^2(x)}Ψ_2(x)+\Psi_1V_1Ψ_2 dx$$
Please notice that $\Psi_I$ relates to region while $\Psi_1$ relates to energy level.
You do have a nontrivial point in that the orthogonality relation for the Laguerre polynomials as they appear in the hydrogenic eigenfunctions,
$$
\int_0^\infty L_{n-\ell-1}^{(2\ell+1)}(2r/n)L_{n'-\ell-1}^{(2\ell+1)}(2r/n') e^{-\frac{n+n'}{nn'}r} r^{2\ell+2} \mathrm dr = 0 \qquad \text{for }n\neq n',
$$
is structurally very different to the standard orthogonality relation,
$$
\int_0^\infty L_{n}^{(2\ell+1)}(r)L_{n'}^{(2\ell+1)}(r) e^{-r} r^{2\ell+1} \mathrm dr = 0 \qquad \text{for }n\neq n',
$$
as provided in e.g. Wikipedia or the DLMF.
As pointed out in the comments, the orthogonality of hydrogenic wavefunctions follows directly from the general Sturm-Liouville theory (they're eigenfunctions of a hermitian Sturm-Liouville operator with different eigenvalues and that's all you should need) so that an elementary proof of the altered orthogonality relation seems a bit of a waste of time.
However, that obviously hasn't stopped other people from looking for it - one suitable example appears to be
A Laguerre Polynomial Orthogonality and the Hydrogen Atom. Charles F. Dunkl. arXiv:math-ph/0011021.
Best Answer
This problem could be done more simply through the application of linear algebra. You want to prove that
$$\langle \psi_1 - \psi_2 | \psi_1 + \psi_2 \rangle = 0$$
The inner product is analogous to the dot product of linear algebra, and it is distributive. Distributing, we find that
$$\begin{aligned} \langle \psi_1 - \psi_2 | \psi_1 + \psi_2 \rangle &= \langle \psi_1 - \psi_2 | \psi_1 \rangle + \langle \psi_1 - \psi_2 | \psi_2 \rangle \\ &= \langle \psi_1 | \psi_1 \rangle - \langle \psi_2 | \psi_1 \rangle + \langle \psi_1 | \psi_2 \rangle - \langle \psi_2 | \psi_2 \rangle \end{aligned} $$
Because $\psi_1$ and $\psi_2$ are orthogonal and normalized, you know $\langle \psi_i | \psi_j \rangle = \delta_{i j}$. Substituting, the above expression evaluates to $1 - 0 + 0 - 1 = 0$, demonstrating that the two vectors are indeed orthogonal.
Your approach - using the integrals - was also valid, and fundamentally similar to mine here. However, by noting that the relation you used ($\langle \psi_1 | \psi_2 \rangle = \int_{-\infty}^{\infty} \! \psi_1^* \psi_2 \, \mathrm{d}x$) satisfied the definition of an inner product, the integrals can be omitted.