How Do Units Affect Hydrogen Wavefunction Normalization in FORTRAN?

[peanutbutter]

I am solving the Hydrogen wavefunction using FORTRAN.
Now using the Euler method, I am given a solution to match which is given by u10(r) = 1.06r*exp(-3.74r) (where unl(r) = rRnl(r) in general) which says it has a normalisation chosen to match what i should get from my code.

Then I use the backwards Euler method getting a more accurate solution with a correct shape and a much larger amplitude. I then normalise this and plot it against the analytic solution given by
u = 1/pi^1/2*(1/ao)^3/2*r*exp(-r/a0) which is apparently also normalised.

I am working in natural units so h-bar = c = 1 and me = 511.7keV the electron mass, which implies that my Bohr radius = 0.26737keV^-1 and I am working in units of r given by keV^-1 also.

so this would imply that the units of u for the analytic solution must be keV^(1/2) however looking at the original equation shouldn't the units be keV^-1 as u = rR(r) where r has units keV^-1 and R(r) is dimensionless, I'm confused. Is my normalisation unit calculated in fortran supposed to have units like the analytic normalisation constant?

On my graphs should I just put no units? :S
Any suggestions would be appreciated, units of wavefunctions has always confused me.

 

[Redbelly98]

Wavefunction units can be tricky, but you can get it from the normalization condition:

∫ Ψ*Ψ dx dy dz = 1​

Since the "1" is dimensionless, and "dx dy dz" is a volume, the wavefunction must have units of volume^-1/2 or length^-3/2