About the standard textbook example for Laplace Equation and Separation of Variables Here is a standard example typically cited to illustrate the technique of separation of variables. Suppose we have a semi-infinite rectangular strip formed by boundaries $x=0, y=0, y=a$ and we are to solve Laplace equation $\Delta f=0$ with boundary conditions (i) $f=0$ when $y=0$ or $y=a$ (ii) $f \to 0$ as $x \to \infty$ (iii) $f\equiv C$ when $x=0$ ($C$ is a constant).
A typical solution by separation of variables will first yield the "solution pieces" $f_n=c_n e^{-n\pi x/a} \sin(n\pi y/a)$ (where $n=1,2,...$) which satisfies (i) and (ii), and then use the method of "handwave" and claim that by linearity, we can consider the candidate 
$f_\infty := \sum_{n=1}^\infty c_n e^{-n\pi x/a} \sin(n\pi y/a)$
and try to find $c_n$ such that (iii) holds.
My question: Yes, the Laplace operator is linear but that only allows for us to piece together finitely many solutions. How do we know that the infinite sum $f_\infty$ is still "differentiable" within the strip? (by "differentiable" here, I meant that the second partial derivatives exist). For example, along the line $x=0$, isn't it the case that $f(0,0)=f(0,a)=0$ but $f(0,y)=C$ for $0<y<a$ so the fitted $f_\infty$ could not be differentiable at $(0,0)$ and $(0,a)$. Are these the only two spots where $f_\infty$ is allowed to be (or turns out to be) non-differentiable?
In addition, this example confuses me because when I learnt the general formulation of a Laplace equation problem, there is the PDE $\Delta f=0$ and boundary conditions, and I always thought that $\Delta f$ should at least make sense on the boundary but in this example, the points $(0,0)$ and $(0,a)$ are doomed to fail. If we completely allow $f$ to be non-differentiable at boundary, then theorems like uniqueness of solution would not hold (e.g. piecewise define $f$ by choosing any $f$ that satisfies $\Delta f=0$ on the interior and separately take the boundary values). So why is it that we consent the poor behavior of the points $(0,0)$ and $(0,a)$ here?
 A: One of the nicest classical results that helps in many cases is Abel's Theorem for the boundary behavior of a power series. In this theorem, one assumes that
$$
             p(z)=\sum_{n=0}^{\infty}c_n z^{n}
$$
converges conditionally or absolutely at some real $z=z_0 > 0$. The conclusion is that $p(x)$ is continuous on $[0,z_0]$ with
$$
             \lim_{z\uparrow z_0} p(x) = p(z_0).
$$
Furthermore, the power series is infinitely differentiable for $|z| < |z_0|$.
To apply this result to your case, note that the following is a power series in $z=e^{-\pi x/a}$:
$$
        f(x,y) = \sum_{n} c_{n}(e^{-\pi x/a})^{n}\sin(n\pi y/a)
$$
Then you know convergence properties when $z=1$ (i.e., where $x=0$) based on properties of the Fourier series $f(0,y)=C=\sum_{n}c_n \sin(n\pi y/a)$. So you get continuity properties of $f(x,y)$ as $x\downarrow 0$, but you won't necessarily get continuity of higher order derivatives unless the derived Fourier series converges as well, which it will normally not do. However, you'll get all orders of derivatives in $x,y$ for $x,y$ strictly inside the open region of the problem.
