I know of no book where this is fully covered for the real case, but I've outlined such an adaption of a common proof given for the complex case. The proof relies on a 1911 result of Herglotz now known as the Herglotz representation theorem for non-negative harmonic functions on the unit disk or half-plane. (See http://en.wikipedia.org/wiki/Herglotz_representation_theorem .)
One book giving a proof using the Herglotz Theorem for the complex case is Gerald Teschl's book Mathematical Methods of Quantum Mechanics. Chapter 3 is The Spectral Theorem, and the required Herglotz Theorem is proved in Appendix A at the end of this chapter, but only in the complex case where $\Im F(z) \ge 0$. (It seems that a full proof for real functions only must be carried out on the unit disk, and mapped to the upper half-plane.) This book is legitimately offered for free for private use. Look for the Download section on Prof. Teschl's web page: http://www.mat.univie.ac.at/~gerald/ftp/book-schroe/index.html . According to the author: "Permission is granted to retrieve and store a single copy for personal use only. If you like this book and want to support the idea of online versions, please consider buying this book."
I'll outline an adapted proof which works for a densely-defined selfadjoint linear operator A on a real Hilbert space $X$. What's interesting is that the proof is constructive, making it particularly interesting for Quantum Mechanics. The proof relies on the Herglotz representation for non-negative harmonic functions on the upper half-plane; Herglotz' representation is fundamental to a large body of holomorphic and harmonic function theory on the unit disk and in the half plane. In it's most general form, the Herglotz representation holds for all non-negative harmonic functions on the unit disk, and, using this representation, one is able to show that all such functions have pointwise a.e. radial limits, and have Poisson integral representations with respect to finite positive Borel measures on the unit circle. The version on the upper half-plane is a little more subtle, but the representation is just as complete.
The idea is to start with a densely-defined selfadjoint linear operator $A : \mathscr{D}(A)\subseteq X\rightarrow X$ on a Hilbert space $X$, and to consider
$$
``\frac{y}{(A-xI)^{2}+y^{2}I}" = y((A-xI)^{2}+y^{2})^{-1}.
$$
This gives rise to a quadratic form $\phi(x,y)$ on $X\times X$:
$$
\phi(x,y)(a,b) = (y((A-xI)^{2}+y^{2}I)^{-1}a,b),\;\;\; x,y\in\mathbb{R},\;y > 0,\;\;\; a,b \in X.
$$
Though not immediately obvious, $(A-xI)^{2}$ is densely-defined, selfadjoint and $(A-xI)^{2}+y^{2}I$ is invertible for all $y \ne 0$.
Moreover, for fixed $a,b\in X$, the above is a non-negative harmonic function of $x, y$ in the open upper half-plane $U$. And, it can be shown that
$$
\lim_{y\uparrow\infty}y\phi(x,y)(a,a) = \|a\|^{2}, \;\;\; a \in X,\\
\phi(x,y)(a,a) \le \|a\|^{2}/y,\;\;\; y > 0.
$$
These are necessary and sufficient conditions in order for there to exist a unique finite positive Borel measure $\mu_{a}$ on $\mathbb{R}$ such that
$$
\phi(x,y)(a,a) = \int_{-\infty}^{+\infty}\frac{y}{(t-x)^{2}+y^{2}}d\mu_{a}(t),\;\;\; \mu_{a}(\mathbb{R})=\|a\|^{2}.
$$
(When $\mu_{a}$ is treated as a non-decreasing bounded function on $\mathbb{R}$, this result is a special case of Herglotz' Theorem.) The uniqueness of $\mu_{a}$ as a measure forces $\mu_{a}$ to inherit the Parallelogram law:
$$
\mu_{a+b}(E)+\mu_{a-b}(E)=2\mu_{a}(E)+2\mu_{b}(E),\;\;\; a,b \in X,
$$
for every Borel subset $E$ of $\mathbb{R}$. And
$$
\mu_{\rho a}(E)= \rho^{2}\mu_{a}(E),\;\;\; \rho \in \mathbb{R},\; a\in X.
$$
Therefore, for a fixed Borel set $E$, $a\rightarrow \mu_{a}(E)$ is a positive bounded quadratic form on $X$. So there is a unique $P(E) \in \mathscr{L}(X)$ such that $P(E)$ is bounded, selfadjoint, positive with $0 \le P(E) \le I$, and satisfies $(P(E)a,a)=\mu_{a}(E)$ for all $a \in X$. Furthermore $P(E)\le P(E')$ if $E\subseteq E'$, and $P(\mathbb{R})=I$.
It follows that
$$
\left(\frac{y}{(A-xI)^{2}+y^{2}I}a,a\right)=\int_{-\infty}^{+\infty}\frac{y}{(t-x)^{2}+y^{2}}d(P(t)a,a),\;\;\; y > 0,\; x\in\mathbb{R},\;\; a \in X.
$$
Basically, this is a Poisson representation of $y/((A-xI)^{2}+y^{2}I)$. One can explicitly integrate (think $\tan^{-1}$) to show that $P$ is determined from $A$ through the integral
$$
\lim_{y\downarrow 0}\frac{1}{\pi}\int_{r}^{s}\left(\frac{y}{(A-xI)^{2}+y^{2}I}a,b\right)dt = \frac{1}{2}(P[r,s]a,b)+\frac{1}{2}(P(r,s)a,b).
$$
Using this, it is possible to show that $P$ is the desired spectral measure for $A$. If $F$ is a compactly supported continuous function on $\mathbb{R}$, then properties of the Poisson integral give an explicit formula for the functional calculus:
$$
F(A)=\lim_{y\downarrow 0}\frac{1}{\pi}\int_{-\infty}^{\infty}f(x)y((A-xI)^{2}+y^{2}I)^{-1}dx.
$$
If $F$ has a few discontinuities, such as when $F=\chi_{[a,b]}$, then the above still exists in the strong operator topology, but not necessarily in the uniform operator topology. Of course the final result is
$$
x \in \mathscr{D}(A) \iff \int_{-\infty}^{\infty}\lambda^{2}d(P(t)x,x) < \infty,
$$
and, in that case,
$$
Ax = \int \lambda dP(\lambda)x,\;\;\; x \in \mathscr{D}(A).
$$
Best Answer
You know that $E(S)x\in\mathcal{D}(A)$ for all $x$ if $S$ is a bounded measurable subset of $\mathbb{R}$. And $$ AE(-R,R]x = \int_{-R}^{R}\lambda dE(\lambda)x. $$ Furthermore, $$ \|AE(-R,R]x\|^2 = \int_{-R}^{R}\lambda^2d\|E(\lambda)x\|^2. $$
If $x\in\mathcal{D}(A)$, then $AE(-R,R]x=E(-R,R]Ax$ converges to $Ax$ as $R\rightarrow\infty$ and, therefore, the above gives $$ x \in\mathcal{D}(A)\implies \int_{-\infty}^{\infty}\lambda^2d\|E(\lambda)x\|^2 < \infty. \tag{*} $$ Conversely, if $\int_{-\infty}^{\infty}\lambda^2 d\|E(\lambda)x\|^2 < \infty$, then $AE(-R,R]x=\int_{-R}^{R}\lambda dE(\lambda)x$ converges as $R\rightarrow\infty$ to some $y$. Because $E(-R,R]x\rightarrow x$ as $R\rightarrow\infty$, and $AE(-R,R]x\rightarrow y$ as $R\rightarrow\infty$, and because $A$ is closed, then $x\in\mathcal{D}(A)$ and $Ax=\int_{-\infty}^{\infty}\lambda dE(\lambda)x$. Hence, the reverse implication is proved for $(*)$, which gives the full characterization of $\mathcal{D}(A)$.