You can download their example from https://web.stanford.edu/~hastie/CASI_files/DATA/kidney.txt and easily replicate their results.
> kidney <- read.table("kidney.txt", header=TRUE)
> str(kidney)
'data.frame': 157 obs. of 2 variables:
$ age: int 18 19 19 20 21 21 21 22 22 22 ...
$ tot: num 2.44 3.86 -1.22 2.3 0.98 -0.5 2.74 -0.12 -1.21 0.99 ...
> fit <- lm(tot ~ age, data=kidney)
> fit$coefficients
(Intercept) age
2.86002680 -0.07858842
As to standard errors, standard errors for fitted values, $\text{se}(\hat{y})$, are differente from coefficients' standard errors, $\text{se}(\hat\beta)$.
The model matrix $X$ is:
> X <- model.matrix(fit)
> head(X)
(Intercept) age
1 1 18
2 1 19
3 1 19
4 1 20
5 1 21
6 1 21
Putting $S=(X^TX)^{-1}$, $\text{cov}(\hat\beta)=\sigma^2_yS$ (see my answer to this question). Given a single fitted value, $\hat{y}_h$ and the corresponding $h$th row of $X$, e.g.
$$y_1=2.44,\qquad x_1=\begin{bmatrix}1 \\ 18\end{bmatrix}$$
the variance of $\hat{y}_h$ is:
$$\text{var}(\hat{y}_h)=\text{var}(x_h^T\hat\beta)=x_h^T\text{cov}(\hat\beta)x_h=x_h^T(S\sigma^2_y)x_h
=\sigma^2_y(x_h^TSx_h)$$
You estimate $\sigma^2_y$ by the residual mean square, RMS, the stardard error of $\hat{y}_h$ is:
$$\text{se}(\hat{y}_h)=\sqrt{RMS(x_h^TSx_h)}$$
and it depends on $x_h$.
When there is only one independent variable,
$$S=(X^TX)^{-1}=\frac{1}{n\sum(x_i-\bar{x})^2}
\begin{bmatrix}\sum x_i^2 & -\sum x_i \\ -\sum x_i & n \end{bmatrix}$$
and
\begin{align*}
x_h^T(X^TX)^{-1}x_h &=\frac{\sum x_i^2-2x_hn\bar{x}+nx_h^2}{n\sum(x_i-\bar{x})^2}=\frac{\sum x_i^2 -n\bar{x}^2+n(x_h-\bar{x})^2}{n\sum(x_i-\bar{x})^2}\\
&=\frac1n+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}
\end{align*}
(Remember that $\sum(x_i-\bar{x})^2=\sum x_i^2-n\bar{x}^2$).
The "extended version of formula (1.2)" (which is just the standard error of a mean) is:
$$\text{se}(\hat{y}_h)=\left[RMS\left(\frac1n+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)\right]^{\frac12}$$
BTW, this is how confidence bands are calculated.
See Kutner, Nachtsheim, Neter & Li, Applied Linear Statistical Models, McGraw-Hill, 2005, §2.4, or Seber & Lee, Linear Regression Analysis, John Wiley & Sons, 2003, §6.1.3.
In R:
> S <- solve(t(X) %*% X)
> RMS <- summary(fit)$sigma^2
> x_h <- matrix(c(1, 20), ncol=1) # first standard error in Table 1.1
> y_h_se <- sqrt(RMS * (t(x_h) %*% S %*% x_h)); y_h_se
[,1]
[1,] 0.2066481
> x_h <- matrix(c(1, 80), ncol=1) # last standard error in Table 1.1
> y_h_se <- sqrt(RMS * (t(x_h) %*% S %*% x_h)); y_h_se
[,1]
[1,] 0.420226
EDIT
If you are interested in the standard error of $\hat{y}_{h(new)}=\hat\alpha+\hat\beta x_{h(new)}$, when $x_{h(new)}$ is a new observation, you do not know what $\hat{y}_h$ would be in a regression on $n+1$ points. Different samples would yield different predictions, so you have to take into account the deviation of $\hat{y}_{h(new)}$ around $\hat{y}_h=\hat\alpha+\hat\beta x_h$:
$$\text{var}[y_{h(new)}-\hat{y}_h]=\text{var}(y_{h(new)})+\text{var}(\hat{y}_h)$$
So the variance of your prediction has two components: the variance of $y$, which you estimate by RMS, and the variance of the sampling distribution of $\hat{y}_h$, $RMS(x_h^TSx_h)$:
$$RMS + RMS\left(\frac1n+\frac{(x_h-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)$$
The "extended version of formula (1.2)" turns into:
$$\text{se}(\hat{y}_{h(new)})=\left[RMS\left(1+\frac1n+\frac{(x_{h(new)}-\bar{x})^2}{\sum(x_i-\bar{x})^2}\right)\right]^{\frac12}$$
See Kutner, Nachtsheim, Neter & Li, Applied Linear Statistical Models, McGraw-Hill, 2005, §2.5, or https://online.stat.psu.edu/stat501/lesson/3/3.3.
