Today I think I can answer the question myself, and considering the renewed interest on it, I shall do so. If anything is found to be mistaken, comments are welcome.
The reason for the question is that Griffiths was imprecise, he gave the idea of the proof, but not the mathematical rigor details. First, let us state better what is the result we want to prove:
Theorem: Let $\mathbf{E}$ be an electrostatic field, i.e. in particular obeying Gauss's law in integral form. Let $S\subset \mathbb{R}^3$ be one charged surface with surface charge $\sigma : S\to \mathbb{R}$ and with normal $\mathbf{n}_0 : S\to T\mathbb{R}^3$. Then if $p\in S$ it holds
$$\lim_{\epsilon\to 0}\mathbf{E}(p+\epsilon \mathbf{n}_0(p))\cdot \mathbf{n}_0(p)-\mathbf{E}(p-\epsilon \mathbf{n}_0(p))\cdot\mathbf{n}_0(p)=\dfrac{\sigma(p)}{\epsilon_0}.$$
We are going to show this using Griffith's idea, but filling in the details. First let's construct the "wafer-thin gaussian pillbox" and use Gauss' law.
Take $p\in S$ a point on the surface. Let $\mathbf{n}_0 : S\to T\mathbb{R}^3$ be the surface's normal vector field. Then let $U\subset S$ be open containing $p$, and consider the following set of points of $\mathbb{R}^3$
$$\mathcal{D}_U(\epsilon)=\{q\in \mathbb{R}^3 : q = q_0 + \lambda \mathbf{n}(q_0),\quad q_0\in U,\lambda \in [-\epsilon,\epsilon]\}$$
In other words: pick a neighborhood of $p$ on the surface. Take every point of the neighborhood and "go out of $S$ from it" in both directions. This is the "gaussian pillbox".
Define also $\mathcal{D}_U^\lambda(\epsilon)$ to be the subset of $\mathcal{D}_U(\epsilon)$ corresponding to one specific $\lambda\in [-\epsilon,\epsilon]$. It should be clear that the area satisfies $A(\mathcal{D}_U^\lambda(\epsilon))=A(U)$.
Now, $\mathcal{D}_U(\epsilon)$ is a three-dimensional volume with a boundary $\partial \mathcal{D}_U(\epsilon)$. Let's analyze its boundary. By construction its boundary can be split up in three parts
$$\partial \mathcal{D}_U(\epsilon)= \mathcal{D}_U^\epsilon(\epsilon)\cup \mathcal{D}_U^{-\epsilon}(\epsilon)\cup \mathcal{W}_U(\epsilon)$$
the first and second terms are respectivelly the top and bottom lids and the third, $\mathcal{W}_U(\epsilon)$ is the wall made up by picking the line bounding $U$ and moving it straight up and down. The wall has area $A(\mathcal{W}_U(\epsilon))=2\ell\epsilon$, where $\ell$ is the lenght of the curve bounding $U$.
Bringing these together, and remembering the only charge inside of $\mathcal{D}_U(\epsilon)$ is that one from the surface, Gauss's law
$$\int_{\partial \mathcal{D}_U(\epsilon)} \mathbf{E}\cdot d\mathbf{a}=\dfrac{1}{\epsilon_0}\int_{U} \sigma dA$$
can be recast as
$$\int_{\partial \mathcal{D}_U^\epsilon(\epsilon)}\mathbf{E}\cdot d\mathbf{a}+\int_{\partial \mathcal{D}_U^{-\epsilon}(\epsilon)} \mathbf{E}\cdot d\mathbf{a} + \int_{\mathcal{W}_U(\epsilon)}\mathbf{E}\cdot d\mathbf{a}=\dfrac{1}{\epsilon_0}\int_U \sigma dA.$$
Now we apply the the mean value theorem for integrals on both sides. This yields $p^\pm \in \mathcal{D}_U^{\pm \epsilon}(\epsilon)$ on the top/bottom lids, $q\in \mathcal{W}_U(\epsilon)$ and $p'\in U$ such that the equation becomes
$$\mathbf{E}(p^+)\cdot \mathbf{n}_0(p^+) A(\mathcal{D}_U^{\epsilon}(\epsilon))-\mathbf{E}(p^-)\cdot \mathbf{n}_0(p^-) A(\mathcal{D}_U^{-\epsilon}(\epsilon))+\mathbf{E}(q)\cdot \mathbf{n}_W(q) A(\mathcal{W}_U^{\epsilon}(\epsilon))=\dfrac{\sigma(p')}{\epsilon_0}A(U).$$
where we have used that the normals of $\mathcal{D}_U^{\pm \epsilon}(\epsilon)$ are just $\mathbf{n}_0$ copied over, and with direction flipped on the bottom lid, and the normal of the wall is $\mathbf{n}_W$. Inserting what we know about the areas gives
$$\mathbf{E}(p^+)\cdot \mathbf{n}_0(p^+) A(U)-\mathbf{E}(p^-)\cdot \mathbf{n}_0(p^-) A(U)+\mathbf{E}(q)\cdot \mathbf{n}_W(q) 2\ell \epsilon =\dfrac{\sigma(p')}{\epsilon_0}A(U).$$
Finally taking $\epsilon\to 0$ yields
$$\lim_{\epsilon \to 0} \mathbf{E}(p^+)\cdot \mathbf{n}_0(p^+) -\mathbf{E}(p^-)\cdot \mathbf{n}_0(p^-) =\dfrac{\sigma(p')}{\epsilon_0}.$$
Now notice that since $p^{\pm} \in \mathcal{D}^{\pm \epsilon}_U(\epsilon)$ it must be $p^+ = p_0 + \epsilon \mathbf{n}_0(p_0)$ and $p^- = p_0' - \epsilon \mathbf{n}_0(p_0')$. This in turn gives a result
$$\lim_{\epsilon \to 0} \mathbf{E}(p_0 + \epsilon \mathbf{n}_0(p_0))\cdot \mathbf{n}_0(p_0) -\mathbf{E}(p_0' - \epsilon \mathbf{n}_0(p_0'))\cdot \mathbf{n}_0(p_0') =\dfrac{\sigma(p')}{\epsilon_0}.$$
Up to now we have left $U$ arbitrary and thus it holds for any such $U$. Consider now $U = B(p,1/n) \cap S$ the open ball centered in $p$ of radius $1/n$. The results will depend on $n$ of course. Thus we shall have $p_0 = x_n$, $p_0' = y_n$ and $p' = z_n$. The equation then becomes
$$\lim_{\epsilon \to 0} \mathbf{E}(x_n + \epsilon \mathbf{n}_0(p_0))\cdot \mathbf{n}_0(x_n) -\mathbf{E}(y_n - \epsilon \mathbf{n}_0(y_n))\cdot \mathbf{n}_0(y_n) =\dfrac{\sigma(z_n)}{\epsilon_0}.$$
Finally take $n\to \infty$ shrinking down to $p$. Since $x_n,y_n,z_n\in B(p,1/n)$ these sequences all converge to $p$. Assuming smoothnes of everything to be able to exchange the limits, we obtain the result
$$\lim_{\epsilon \to 0} \mathbf{E}(p + \epsilon \mathbf{n}_0(p))\cdot \mathbf{n}_0(p) -\mathbf{E}(p- \epsilon \mathbf{n}_0(p))\cdot \mathbf{n}_0(p) =\dfrac{\sigma(p)}{\epsilon_0},$$
as we wished to prove.