Fourier transform of a wavefunction After solving the Schrodinger equation numerically by finite difference method for some potential.
I can plot the wavefunction with the distance from nucleus and I get the right theoretical shape. 
Now if I take the Fourier transform of the wavefunction I get some numbers but I do not have the corresponding k values to plot the momentum space wavefunction. Where is my problem? 
Again, I have $r$ and $\psi(r)$, I take FFT of $\psi(r)$ and I apparently get $\psi(k)$, if I want to plot $\psi(k)$ against $k$, where do I get $k$ from? 
Did I miss something elementary?
 A: Analytics
First note that if you are just naively using some black-box FFT function, you are likely not computing what you think you are computing. You have a function that depends only on the length of a 3-d position, $\psi(r)$. Its Fourier transform in 3-d can be written in spherical coordinates as
$$ \psi(k) = \int \psi(r) e^{-i \vec r\cdot \vec k} d^3 r = \int_0^\infty\int_0^\pi\int_0^{2\pi} \psi(r) e^{-i k r \cos\theta} r^2 \sin\theta d\phi d\theta dr
$$
where I've defined $\theta$ to be the angle between $\vec k$ and $\vec r$. You can perform the angular integrals analytically. Doing so, you will get
$$ \psi(k) = \frac{4\pi}{k}\int_0^\infty \psi(r) r \sin(kr) dr $$
Evidently, the Fourier transform in 3-d of $\psi(r)$ can instead be obtained from the sine transform of $r\psi(r)$. However, this is only one-sided ($0$ to $\infty$). Most Fourier transform routines are implemented based on two-sided Fourier transforms, where the variables run from $-\infty$ to $\infty$. If we assume $\psi(r)$ is an even function, $\psi(r)=\psi(-r)$, then we can say that 
$$
\psi(k) = \frac{2\pi}{k}\int_{-\infty}^\infty r\psi(r)\sin(kr)dr 
$$
Some FFT libraries include discrete sine transforms. If yours does not, you can use the fact that $\sin(x) = -\mathrm{Im}[e^{-ix}]$ to express your transform as 
$$ \psi(k) = -\frac{2\pi}{k} \mathrm{Im} \int_{-\infty}^\infty r\psi(r) e^{-ikr} dr$$
This now involves an ordinary (complex-valued) Fourier transform. A key point is that we've reformulated the transformation between the spherical variables $r$ and $k$ in such a way that it "looks like" a transformation between a Cartesian coordinate $x$ and wavenumber $k_x$.
Numerics
The spacing between your tabulated $r$-values sets the largest $k$ resolved by the Fourier transform. Likewise, the largest $r$-value sets $k$ spacing you can resolve. If it were a 1-d transformation $x\to k_x$, these would be related by
$$ \Delta x \cdot k_{\max} = 2\pi \\ x_{\max} \cdot \Delta k = 2\pi$$
Or, more generally, 
$$ \Delta x \cdot \Delta k = \frac{2\pi}{N}$$
where $N$ is the number of tabulated values. This assumes that $x$ runs from $-x_{\max}$ to $x_{\max}$. Your case is different because your $r$-values only run from $0$ to $r_{\max}$ in $N$ points. So when you extend to negative $r$ as done above, it is as if you double the number of points to $2N$. So your resolveable wavenumbers are given by
$$\Delta r \cdot \Delta  k = \pi/N \implies \Delta k = \pi/r_{\max}$$
So your $k$ values are linearly spaced with the above $\Delta k$.
