Laplacian of 1/r satisfies laplacian equation? We know that 
$$ \nabla^2 \left( \frac{1}{r} \right) = -4 \pi \delta(r) \tag{1}$$
and that the general solution to the laplace equation $\nabla^2\Psi = 0$
may be expanded as
$$\Psi=\sum_l \left(A_lr^l+B_lr^{-(l+1)}\right)P_l(\cos\theta) \tag{2}$$
where I've kept it simple by assuming azimuthal symmetry. 
If we take the $l=0$ term and apply the laplacian operator, we get
$$\nabla^2\Psi_{l=0} = -4\pi B_0\delta(r) \, , \tag{3}$$
so the laplacian is only zero for $r \ne 0$. However, I thought the general solution $(2)$ of the Laplace equation was supposed to hold for all $r$ and $\theta$. Yet $r=0$ it is clearly not a solution. What gives? 
Granted I can understand that a $\frac{1}{r}$ potential is equivalent to that of a point charge, in which the Laplace equation would no longer apply, but I thought the general solution would hold for all regions of interest in which we know there is no charge. 
 A: Equation $(2)$ is indeed a general solution, but that doesn't mean that all the $A_l$ and $B_l$ have to be nonzero all the time.
For a problem in which $r=0$ is part of the domain the $B_l$ coefficients are zero, else the potential diverges at $r=0$ due to the $r^{-(l+1)}$ functions.
In the case of the point charge, the point $r=0$ is not part of the problem's domain (this is a little detail we normally sort of ignore) and we can have $B_l$ terms.
As you know, for the point charge there's exactly one, $B_1$.
Also for the point charge, the $A_l$ are all zero, as explained below.
Consider a problem with a spherical region cut out of a conducting block, as shown in Figure 1.
The boundary conditions here are that the potential is $V_\text{inner}$ on the inner conductor and $V_\text{outer}$ on the outer one.
The potential in the gap area can contain contributions from the $A_l$ terms and the $B_l$ terms because neither of those diverge in that area.

Figure 1
However, in a situation where the potentials extend to $r=\infty$, as shown in Figure 2, there's a boundary condition that the potential at infinity is zero.$^{[a]}$
Therefore, the $A_l$ terms must all be zero, otherwise the boundary conditions are not met due to the $r^l$ functions.

Figure 2
$[a]$: Putting the potential to zero at infinity is a convenient choice. You can pick any other constant but that doesn't modify the physics.
A: If you solve for just inside a sphere you might need to throw out your B terms.
If you solve for just outside a sphere you might need to throw out your A terms.
But if you are solving for the region between two spherical shells you might need both.
That's why it is general, because in general you need them.
A: The standard math answer is that the delta function fails certain regularity conditions in order for the equation you wrote to make sense. 
However, it is more meaningful to step back and think about where the Laplace equation comes from physically. Specifically, it arises when a quantity (molecular concentration, electric field, thermal energy, whatever) is diffusing in a manner proportional to it's gradient, while also being subject to a conservation law.
Thus the natural statement of the Poisson equation is not:
$$\Delta \Psi = f,$$
but rather, for any open region $\omega$ no matter how small, the outflow minus the inflow of $\Psi$ through the boundary equals the difference between the sources and sinks inside. Symbolically,
$$\underbrace{-\int_{\partial\omega}\nabla \Psi \cdot n dS}_{\text{outflow}-\text{inflow}} = \underbrace{\int_\omega f d\bf{x}}_{\text{sources}-\text{sinks}}, \quad \text{ for all open sets }\omega. $$
Then one uses the divergence theorem, takes the size of $\omega$ to zero, and makes some regularity assumptions on $\Psi$ and $f$ to arrive at the standard equation.
But look at the physical (inflow - outflow) condition - this is a condition based on integrating over a small region. You can verify this condition by computing integrals of delta functions over small regions, rather than "evaluating" them at points. This is a condition that your function does satisfy.
