This is the same idea as timur's answer but with more details and less generality. A frequent test problem in numerical analysis is the Poisson equation $-\Delta u = 1$ on the L-shaped domain
$\Omega = ([-1,1] \times [-1,1]) \setminus ([-1,0] \times [-1,0])$
with homogeneous Dirichlet boundary conditions: $u = 0$ on $\partial\Omega$. The solution has a singularity at the origin: it is continuous but not differentiable. More precisely, close to the origin we have
$u(r,\theta) \approx r^{2/3} \sin \frac{2\theta+\pi}{3}$
in polar coordinates, according to equation (1.6) in http://eprints.ma.man.ac.uk/894/02/covered/MIMS_ep2007_156_Sample_Chapter.pdf (sample chapter from Elman, Silvester and Wathen, Finite Elements and Fast Iterative Solvers, Oxford University Press, 2005).
Added: I don't know the details and I don't have time to do the necessary computations, but I think that you can solve the PDE by converting the Laplacian to polar coordinates and applying separation of variables. I imagine that you get that
$u(r,\theta) = r^{2n/3} \sin \frac{2n}{3} (\theta + \frac{1}{2}\pi)$
with $n$ a positive integer satisfies the boundary conditions at $r=0$ and $\theta=-\pi/2$ and $\theta=\pi$ (as Dorian comments below, these are all harmonic functions, so there must be something else). Then take a linear combination of those to match the conditions on the rest of the boundary of the L-shaped domain. Close to the origin, the $n=1$ term dominates. Perhaps somebody else can confirm / amend?