Consider the finite 1D wedge-shaped potential well given by
$$V(x)=V_0\left(\frac{|x|}{a}-1\right) \hspace{10pt}\mathrm{for}\hspace{3pt} |x|<a;\hspace{6pt}V(x)=0 \hspace{10pt}\mathrm{for}\hspace{3pt} |x|>a.$$
I'm trying to find the reflection and transmission coefficients for a stream of electrons coming from $x=-\infty$ (so $E>0$, no bound states involved). In order to do that, I've divided the domain into 4 parts (like in the picture), solved the Schrödinger equation on each one, got a linear system for the coefficients of each wavefunction and solved it. Meanwhile I've also solved the problem numerically, and got pretty expected results, like the one below.
All was well untill I plugged the analytic solution for the transmission coefficient into Python to render a graph and got this. $E$ is in eV, $x$ and $a$ in nm.
Now, I've gone over my calculations a dozen times, and I still can't spot a mistake, so I want to know wheather it's possible to have a transmission coefficient greater than one for certain values of energy, like in this scenario.
I'm aware that this breaks conservation of energy (red alert!), and I'm still hoping that this is due to a numerical error, but I want to hear some other ideas regarding this phenomenon. I've found this as a related question, but I'm not really satisfied with the answers. Why doesn't this phenomenon happen in a rectangular well? Should I discard the $T>1$ solutions and leave energy gaps in the graph? That seems pretty arbitrary arbitrary to me, but it may be the case here.
EDIT: The comments suggest to give an insight into my analytical solution.
Some initial substitutions: $k=\sqrt{\frac{2mE}{\hbar^2}}$, $k_0=\sqrt{\frac{2mV_0}{\hbar^2}}$, $q_0=(k_0a)^{2/3}$, $\epsilon^2=\frac{E}{V_0}$.
The time-independant Schrödinger equation for the wavefunction on each section reads
$$\frac{d^2\psi_I}{dx^2}+k^2\psi_I=0$$
$$\frac{d^2\psi_{II}}{dx^2}+(-\frac{x}{a}+\epsilon^2-1)k_0^2\psi_{II}=0$$
$$\frac{d^2\psi_{III}}{dx^2}+(\frac{x}{a}+\epsilon^2-1)k_0^2\psi_{III}=0$$
$$\frac{d^2\psi_{IV}}{dx^2}+k^2\psi_{IV}=0$$
Solutions to $I$ and $IV$ may be immediatly given as
$$\psi_I(x)=Ae^{ikx}+Be^{-ikx}$$
$$\psi_{IV}(x)=Ge^{ikx},$$
the first corresponding to the incoming and reflected wave, and the other to the transmitted wave.
Introducing the change of variables
$$\xi(x)=q_0(\frac{|x|}{a}+\epsilon^2-1)$$
transforms the other two equations into Airy differential equations:
$$\frac{d^2\psi_{II}}{d\xi^2}=-\xi\psi_{II}$$
$$\frac{d^2\psi_{III}}{d\xi^2}=-\xi\psi_{III}$$
with general solutions in the form
$$\psi_{II}(x)=C\mathrm{Ai}(-\xi)+D\mathrm{Bi}(-\xi)$$
$$\psi_{III}(x)=E\mathrm{Ai}(-\xi)+F\mathrm{Bi}(-\xi),$$
where $\mathrm{Ai}(z)$ and $\mathrm{Bi}(z)$ are the standard Airy functions of the first and second kind.
Using the fact that both $\psi(x)$ and $\psi'(x)$ are continuous, this gives the linear system:
$$\left(\begin{matrix}
e^{-ika} & e^{ika} & -\mathrm{Ai}(-\xi_a) & -\mathrm{Bi}(-\xi_a) & 0 & 0 \\
-\frac{ika}{q_0}e^{-ika} & \frac{ika}{q_0}e^{ika} & \mathrm{Ai}'(-\xi_a) & \mathrm{Bi}'(-\xi_a) & 0 & 0 \\
0 & 0 & \mathrm{Ai}(-\xi_0) & \mathrm{Bi}(-\xi_0) & -\mathrm{Ai}(-\xi_0) & -\mathrm{Bi}(-\xi_0) \\
0 & 0 & \mathrm{Ai}'(-\xi_0) & \mathrm{Bi}'(-\xi_0) & \mathrm{Ai}'(-\xi_0) & \mathrm{Bi}'(-\xi_0) \\
0 & 0 & 0 & 0 & \mathrm{Ai}(-\xi_a) & \mathrm{Bi}(-\xi_a) \\
0 & 0 & 0 & 0 & \mathrm{Ai}'(-\xi_a) & \mathrm{Bi}(-\xi_a)
\end{matrix} \right)\left(\begin{matrix}
A \\
B \\
C \\
D \\
E \\
F
\end{matrix} \right)=\left( \begin{matrix}
0 \\
0 \\
0 \\
0 \\
Ge^{ika} \\
-G\frac{ika}{q_0}e^{ika}
\end{matrix} \right) $$
where $\xi_a\equiv\xi(x=a)=\xi(-a)=q_0\epsilon^2$ and $\xi_0\equiv\xi(0)=q_0(\epsilon^2-1)$ appear as shorthand notations. Prime corresponds to $\frac{d}{dx}$. The currents are given by
$$J_I=\frac{\hbar}{m}\Im\left(\psi_I^*\frac{d\psi_I}{dx}\right)=\frac{\hbar}{m}(|A|^2-|B|^2)=J_{in}-J_{ref}$$
$$J_{IV}=\frac{\hbar}{m}\Im\left(\psi_{IV}^*\frac{d\psi_{IV}}{dx}\right)=\frac{\hbar}{m}|G|^2=J_{tr}$$
Since $T+R=1$, it follows that $|A|^2-|B|^2=|G|^2$. $G$ is actually undetermined and is used as a free coefficient (every other coefficient can be expressed as $L=G\cdot blabla$), so it may be freely set to be 1. It follows that $T=\frac{J_{tr}}{J_{in}}=\frac{1}{|A|^2}$ and $R=1-\frac{1}{|A|^2}$. Hence, only $A$ is needed to calculate the transmission and reflexion coefficients. It may be calculated from the above system using Cramer's rule.
My result:
\begin{equation*}\label{parametara}
\begin{split}
A&=\frac{iq_0\pi^2e^{2ika}}{ka} \left[\left(\left(\frac{ka}{q_0}\right)^2\mathrm{Ai}^2(-\xi_a)+\mathrm{Ai}'^2(-\xi_a)\right)\mathrm{Bi}(-\xi_0)\mathrm{Bi}'(-\xi_0)\right.\\ &+\left(\left(\frac{ka}{q_0}\right)^2\mathrm{Bi}^2(-\xi_a)+\mathrm{Bi}'^2(-\xi_a)\right)\mathrm{Ai}(-\xi_0)\mathrm{Ai}'(-\xi_0)\\
&\left.-\left(\left(\frac{ka}{q_0}\right)^2\mathrm{Ai}(-\xi_a)\mathrm{Bi}(-\xi_a)+\mathrm{Ai}'(-\xi_a)\mathrm{Bi}'(-\xi_a)\right)(\mathrm{Ai}(-\xi_0)\mathrm{Bi}'(-\xi_0)+\mathrm{Ai}'(-\xi_0)\mathrm{Bi}(-\xi_0)) \right],
\end{split}
\end{equation*}
where use was made of the Wronskian for Airy functions:
$$\mathcal{W}[\mathrm{Ai}(z),\mathrm{Bi}(z)]=\pi^{-1}.$$
The next step was to plug the absolute square of $A$ into expressions for $T$ and $R$, which rendered the problematic graph.
Best Answer
I haven't studied your solution in detail, but, for what it's worth, the formulas for the solutions of Airy equation differ depending on the sign of the coefficient. Have you checked that your solutions are applicable independent of the value of $E$? You may wish to check if your solutions give correct results for $V_0=0$.