On page 97 in the book, how did Elias inferred that instead of the factor $(1+2^k|y|/t)^N$ we must instead insert the factor: $\frac{t^L (\epsilon+ 2^{-k}t+\epsilon |x-y|)^L}{(\epsilon+t+\epsilon|x|)^{-L}(2^{-k}t)^{-L}}(1+2^k|y|/t)^N$?; I think that in the denominator the $-L$ should be plainly just $L$.
https://books.google.co.il/books?id=ljcOSMK7t0EC&printsec=frontcover#v=onepage&q&f=false
It's at the bottom of page 97.
I believe the idea of Stein's argument is correct--it is Stein, after all--but I'm having trouble completing the details using the suggestion as written. If you'll allow me to change these around a bit, I believe we can arrive at the desired conlusion. Recall that we want to show that $$\left\|\sup_{\Psi\in S_{\mathcal{F}}}M_{\Psi}^{\epsilon,L}(f)\right\|_{L^{p}}\leq c_{L}\left\|M_{\Phi}^{\epsilon,L}(f)\right\|_{L^{p}},$$ where the bound $c_{L}$ is independent of $\epsilon$ and may depend on $L$. To prove this we "mimic" the argument used to prove (12) on pg. 95.
Recall Lemma 2 on pg. 93. For $t>0$, we have by the lemma that $$\left|(f\ast\Psi_{t})(y)\right|\leq\sum_{k=0}^{\infty}\left|(f\ast\Phi_{2^{-k}t}\ast\eta_{t}^{(k)})(y)\right|$$ Fix $x\in\mathbb{R}^{n}$ and let $y\in\mathbb{R}^{n}$ be such that $\left|x-y\right|<t<\epsilon^{-1}$, where $\epsilon>0$ is a parameter as in Stein. Then \begin{align*} \left|(f\ast\Psi_{t})(y)\right|\dfrac{t^{L}}{(\epsilon+t+\epsilon\left|y\right|)^{L}}&\leq\sum_{k=0}^{\infty}\underbrace{\int_{\mathbb{R}^{n}}\left|(f\ast\Phi_{2^{-k}t})(y-z)\right||\dfrac{t^{L}}{(\epsilon+t+\epsilon\left|y\right|)^{L}}t^{-n}\left|\eta^{(k)}(\frac{z}{t})\right|dz}_{:=I_{k}}\\ \end{align*} Multiply the RHS by \begin{align*} 1&=\dfrac{(2^{-k}t)^{L}}{(\epsilon+2^{-k}t+\epsilon\left|(y-z)\right|)^{L}}\left(\dfrac{(2^{-k}t)^{L}}{(\epsilon+2^{-k}t+\epsilon\left|(y-z)\right|)^{L}}\right)^{-1}\\ &\cdot\left(1+\dfrac{\left|x-(y-z)\right|}{2^{-k}t}\right)^{-N}\left(1+\dfrac{\left|x-(y-z)\right|}{2^{-k}t}\right)^{N} \end{align*}
Since we are going to eventually let $\epsilon\downarrow 0$, we may assume that $\epsilon<1$. We define an $\epsilon$-variant of the tangential maximal function $M_{N}^{**}$ associated to $\Phi$ by \begin{align*} M_{N}^{\epsilon,**}f(x)&:=\sup_{{y'\in\mathbb{R}^{n}}\atop{0<t<\epsilon^{-1}}}\left|(f\ast\Phi_{t})(x-y')\right|\dfrac{t^{L}}{(\epsilon+t+\epsilon\left|x-y'\right|)^{L}}\left(1+\dfrac{\left|y'\right|}{t}\right)^{-N}\\ &=\sup_{{y'\in\mathbb{R}^{n}}\atop{0<t<\epsilon^{-1}}}\left|(f\ast\Phi_{t})(y')\right|\dfrac{t^{L}}{(\epsilon+t+\epsilon\left|y'\right|)^{L}}\left(1+\dfrac{\left|x-y'\right|}{t}\right)^{-N} \end{align*}
We have the simple estimate \begin{align*} \left(1+\dfrac{\left|z\right|}{t}\right)\left(\epsilon+t\epsilon+\left|y\right|\right)&=\epsilon+t+\epsilon\left|y\right|+\dfrac{\epsilon\left|z\right|}{t}+\left|z\right|+\dfrac{\epsilon\left|z\right|\left|y\right|}{t}\\ &\geq\epsilon+2^{-k}t+\epsilon\left|y\right|+\epsilon\left|z\right|\\ &\geq\epsilon+2^{-k}t+\epsilon\left|y-z\right| \end{align*} where the ultimate inequality is just triangle inequality. Together with the fact that $\left|x-y\right|<t$, we obtain the estimate \begin{align*} I_{k}&\leq t^{-n}\int_{\mathbb{R}^{n}}(M_{N}^{\epsilon,**}f)(x)\left(1+\dfrac{2^{k}\left|x-(y-z)\right|}{t}\right)^{N}\dfrac{t^{L}(\epsilon+2^{-k}t+\epsilon\left|(y-z)\right|)^{L}}{(2^{-k}t)^{L}(\epsilon+t+\epsilon\left|y\right|)^{L}}\left|\eta^{(k)}(\frac{z}{t})\right|dz\\ &\leq 2^{kL}t^{-n}(M_{N}^{\epsilon,**}f)(x)\int_{\mathbb{R}^{n}}\left(1+2^{k}+\frac{2^{k}\left|z\right|}{t}\right)^{N}\left(1+\dfrac{\left|z\right|}{t}\right)^{L}\left|\eta^{(k)}(\frac{z}{t})\right|dz\\ &\lesssim_{N} 2^{k(L+N)}(M_{N}^{\epsilon,**}f)(x)\int_{\mathbb{R}^{n}}\left(1+{\left|z\right|}\right)^{N+L}\left|\eta^{(k)}(z)\right|dz \end{align*} where the last inequality is just dilation invariance and pulling out the $2^{k}$ from the first factor of the integrand. If $\left\|\eta^{(k)}\right\|_{\alpha,\beta}\lesssim_{L}2^{-k(N+L+1)}$, then $$I_{k}\lesssim_{L}2^{-k}(M_{N}^{\epsilon,**}f)(x)$$ As Stein notes, this holds if $\Psi$ belongs to a suitably chosen $\mathcal{S}_{\mathcal{F}}$. Therefore for all $\Psi\in\mathcal{S}_{\mathcal{F}}$, we have the estimate $$\left|(f\ast\Psi_{t})(y)\right|\lesssim_{L}\sum_{k=0}^{\infty}2^{-k}(M_{N}^{\epsilon,**}f)(x)\leq (M_{N}^{\epsilon,**}f)(x),\quad\forall 0<t<\epsilon^{-1}, \forall \left|x-y\right|<t$$ Taking the supremum over all such $y$, we conclude that $$(\mathcal{M}_{\Psi}^{\epsilon,L}f)(x)\lesssim_{L} (M_{N}^{\epsilon,**}f)(x),\quad\forall\Psi\in\mathcal{S}_{\mathcal{F}},$$ which implies that $$\sup_{\Psi\in\mathcal{S}_{\mathcal{F}}}(M_{\Psi}^{\epsilon,L}f)(x)\lesssim_{L} (M_{N}^{\epsilon,**}f)(x)$$ To complete the proof, we follow the proof of Lemma 1 on pg. 92. replacing $F_{a}^{*}$ and $F^{*}$ by $\epsilon$-variants. Let me know if you would like details on this point. Taking $N>n/p$, we conclude that $$\left\|\sup_{\Psi\in\mathcal{S}_{\mathcal{F}}}M_{\Psi}^{\epsilon,L}f\right\|_{L^{p}}\lesssim_{L}\left\|M_{N}^{\epsilon,**}f\right\|_{L^{p}}\lesssim_{N ,p,L}\left\|M_{\Phi}^{\epsilon,L}f\right\|_{L^{p}}$$