I will post it as an answer because it's a bit long, although just technical.

We're looking at the Taylor series -
$$ \sum_{n} \frac{1}{n!}E^{Y}_{n}(z,s_0)(s-s_0)^{n} $$
notice that the capital $Y$ won't cause any trouble, we can assume that we are in the compact part away from the cusp.

The $L^2$ estimate proved in Kubota shows that the partial sums of this sequence converges weakly to this sum (which is in $L^2$ due to the truncation).
Now it is enough to show that the partial sums are (generalized-)eigenfunctions and use thm 1.3.4 in Kubota.

So we just apply the Laplacian to it.
$$\Delta \sum_{n} \frac{1}{n!}E^{Y}_{n}(z,s_0)(s-s_0)^{n} = \sum_{n}\frac{1}{n!}[\frac{\partial^{n}}{\partial s^{n}}\Delta E(z,s)]_{s=s_0}(s-s_0)^{n} $$
Because the differentiation is over $s$ and the Laplacian cares only about $x,y$ and this is some instance of Schwartz' lemma of multivariable calculus.

Now $\Delta E(z,s)=s(s-1)E(z,s)$, (I'm using Kubota normalization, where $\Delta$ is positive and not negative).

Therefore we'll calculate $\frac{\partial^{n}}{\partial s^{n}} s(s-1)E(z,s)$.
We'll calculate it by generalized Leibnitz rule.

$$\frac{\partial^{n}}{\partial s^{n}} s(s-1)E(z,s)=\sum_{k=0}^{n}\binom{n}{k}[s(s-1)]^{(k)}E(z,s)^{(n-k)}.$$

Notice that the first function is polynomial in $s$ of deg. $2$, hence we have that
$$\sum_{k=0}^{n}\binom{n}{k}[s(s-1)]^{(k)}E(z,s)^{(n-k)}=\sum_{k=0}^{2}\binom{n}{k}[s(s-1)]^{(k)}E(z,s)^{(n-k)}.$$

By explicitly computing $k=0,1,2$ we have -
$$\sum_{k=0}^{2}\binom{n}{k}[s(s-1)]^{(k)}E(z,s)^{(n-k)}=s(s-1)E(z,s)^{(n)}+n(2s-1)E(z,s)^{(n-1)}+n(n-1)E(z,s)^{(n-2)}.$$

This is of-course not applicable to the case where $n=0,1$, which separate explicit calculation is needed.

Now by grouping the proper derivatives of the Eisenstein series together you can see that
the functions are indeed eigenfunctions of the Laplacian.
$$\frac{s(s-1)}{n!}E(z,s)^{n}+\frac{(n+1)(2s-1)}{(n+1)!}E(z,s)^{(n+1)-1}+\frac{(n+2)(n+1)}{(n+2)!}E(z,s)^{(n+2)-2}=\frac{1}{n!}(s(s-1)+(2s-1)+1)E(z,s)^{(n)}$$

This doesn't apply to the top $2$ powers (you won't sum over the three terms as done above), and I can't seem to fix it in any reasonable way.

Nevertheless, you might still, somehow, use the machinery of thm 1.3.4.
The limit, as the computation shows, is an ($L^2$-) eigenfunction of the Laplacian.

Now one needs to show that this convergence happens in a certain Sobolov space, and not just in $L^2(\Gamma \backslash H)$. If so, then the dirac measure defines a bounded functional on the Sobolov space (by say convolution against the smooth vectors, which you will always have in hand by a theorem due to Harish-Chandra, see Borel's book Ch.2) and then you will get pointwise convergence out of weak convergence.

Now, unless I'm completely wrong in the above calculations, I believe that the original proof in Kubota (and which might be also the one in Jacquet-Gelbart) was flawed, because of this iconic error that Paul is mentioning in his notes, the trucnted Eisenstein series is not an eigenfunction of the Laplacian - http://www.math.umn.edu/~garrett/m/v/iconic_error.pdf

1more comment