Determining power spectral density from a sample
In Matlab, you are dealing with sampled signals, so your frequency variable goes from $0$ to Nyquist. This is why the x-axis is labeled in normalized frequency that goes from $0$ to $\pi$. If you like you can use [P, f] = pwelch(___, fs) where fs is your sampling rate, and then the function will return the power spectral density (PSD) estimate P as a function of analog frequencies in the vector f. See example here.

pwelch is not the only method for estimating PSD. The algorithm choice will depend on your application. A nice summary of other competing algorithms can be found in the answer to this question.

As for these integrals, the limits should go from $-\infty$ to $\infty$ for continuous signals. But if you know that the signal is bandlimited to some range of frequencies $(f_1, f_2)$ then you can get away with integrating from $f_1$ to $f_2$. Recall the definition of PSD is inspired by continuous time signals $x(t)$ with finite energy:
$$
E = \int_{-\infty}^{\infty} |x(t)|^2 dt &lt; \infty
$$
and denoting $X(f)$ to be the Fourier transform of $x(t)$, by Parseval's theorem we have:
$$
\int_{-\infty}^{\infty} |x(t)|^2 dt = \int_{-\infty}^{\infty} |X(f)|^2 df
$$
This leads to the intuition that $|X(f)|^2$ is the "density" of "spectral energy" around the frequency $f.$

Requiring that the PSD integrate to 1 implies that your signals of interest are unit energy. For a non-unit-energy finite energy signal, you can simply divide the signal $x(t)$ by $\sqrt{E}$ to transform it into a unit energy signal.

For implementation purposes, all your signals and transforms are discrete. So you will have to approximate these integrals with sums.