Wavefunctions are found by solving the time-independent Schrödinger equation, which is simply an eigenvalue problem for a well-behaved operator:
$$ \hat{H} \psi = E \psi. $$
As such, we expect the solutions to be determined only up to scaling. Clearly if $\psi_n$ is a solution with eigenvalue $E_n$, then
$$ \hat{H} (A \psi_n) = A \hat{H} \psi_n = A E_n \psi_n = E_n (A \psi_n) $$
for any constant $A$, so $\psi$ can always be rescaled. In this sense, there is no physical meaning associated with $A$.
To actually choose a value, for instance for plotting, you need some sort of normalization scheme. For square-integrable functions, we often enforce
$$ \int \psi^* \psi = 1 $$
in order to bring the wavefunction more in line with the traditional definition of probability (which says the sum of probabilities is $1$, also an arbitrary constant).
In your case,
$$ \psi(x) =
\begin{cases}
\psi_\mathrm{I}(x), & x < -\frac{d}{2} \\
\psi_\mathrm{II}(x), & \lvert x \rvert < \frac{d}{2} \\
\psi_\mathrm{III}(x), & x > \frac{d}{2}.
\end{cases} $$
Thus choose $A$ such that
$$ \int_{-\infty}^\infty \psi(x)^* \psi(x) \ \mathrm{d}x \equiv \int_{-\infty}^{-d/2} \psi_\mathrm{I}(x)^* \psi_\mathrm{I}(x) \ \mathrm{d}x + \int_{-d/2}^{d/2} \psi_\mathrm{II}(x)^* \psi_\mathrm{II}(x) \ \mathrm{d}x + \int_{d/2}^\infty \psi_\mathrm{III}(x)^* \psi_\mathrm{III}(x) \ \mathrm{d}x $$
is unity.
If you happen to be in the regime $E > W_p$, then $\mathcal{K}$ will be imaginary, $\psi_\mathrm{I}$ and $\psi_\mathrm{III}$ will be oscillatory rather than decaying, and the first and third of those integrals will not converge. You could pick an $A$ that conforms to some sort of "normalizing to a delta function," but there are many different variations on this, especially for a split-up domain like this. In that case I would recommend picking an $A$, if you really have to do it, based on some other criterion, such as $\max(\lvert \psi \rvert) = 1$ or something.
Personally, I’d do the approximation in the exponential,
$$e^\epsilon \approx 1 + \epsilon + \frac{\epsilon^2}{2!} + \cdots, $$
before doing the integral.
\begin{align}
P &= \frac4{a_0^3} \int_0^R \mathrm dr\, r^2 \exp\frac{-2r}{a_0}
\\
\frac{a_0^3}{4} P
&\approx
\int_0^R\mathrm dr\,r^2 - \frac2{a_0}\int_0^R \mathrm dr\,r^3
\\&= \frac{R^3}{3} - \frac2{a_0}\frac{R^4}{4}
\\&= \frac{R^3}{3} \left(
1 - \frac32 \frac R{a_0}
\right)
\\
P &\approx \frac43 \left(\frac{R}{a_0}\right)^3
\end{align}
Throwing away the second term basically says “this volume is small enough we can assume the exponential is secretly a constant.”
This result suggests you are getting zeros because you’re only keeping terms that are quadratic in $R/a_0$, but apparently they all get cancelled out.
Best Answer
The plots you see in the Wiki images, are, as their title suggests, probability density plots ($\psi^2$, for Real $\psi$). Light shaded areas represent high probability areas, darker areas lower probability density.
Furthermore, they've been sliced, e.g. with an $x,y$-plane. For radially symmetric functions like $\psi_{2,0,0}$ (aka the $2s$ orbital) the specific slicing plane doesn't even matter: the $\psi^2$ value is the same in all directions as it only depends on $r$.
There's little to be gained from these plots for radially symmetric wave functions with regards to the simpler radial distribution plots like these.
If you really do want to plot radially non-symmetric probability density functions (e.g. for $2p$, 3$d$ or $4f$ orbitals) you'll need to find a tool to generate contour plots. I imagine advanced functions of Mathlab allow to do that.
Here's a *contour plot* of $\psi^2_{200}$ for $z=0$ (the $x,y$-plane):
The white areas are those of high probability density. The plot clearly shows the 'shell-like' structure of the hydrogen $2s$ orbital.
The plot was obtained by means of Wolfram alpha, by plotting the square of:
$$\psi=\frac{1}{4\sqrt{2\pi}a_0 ^\tfrac{3}{2}} \left(2-\frac{\sqrt{x^2+y^2}}{a_0}\right)e^{-\tfrac{\sqrt{x^2+y^2}}{2a_0}}$$
For $a_0=1$.
Here's a excellent resource for visualised wave functions: the Orbitron