To supplement the existing answers, I want to mention a useful application of elliptic regularity for separation of variables in particular.
Let $L$ be a second order elliptic linear operator and consider the PDE,
$$ Lu =f\ \text{ in }\ \Omega, \quad u = 0\ \text{ on }\ \partial\Omega. $$
Here $\Omega$ is a sufficiently nice space parametrised by coordinates $r,x,y$ (for simplicity $\dim \Omega =3$). Assuming the forms are explicitly known, we may try and solve this by separating variables. Doing so, we obtain a solution of the form,
$$ u(x,y) = \sum_{i,j=1}^{\infty} a_{ij}(r)\phi_i(x)\psi_j(y). $$
where $\phi_i,\psi_j$ are an $L^2$-orthonormal basis of eigenfunctions for $L$, satisfying the boundary conditions.
Now assuming the coefficients $a_{ij}$ decay sufficiently quickly, one can show that $u \in H^1(\Omega)$ with weak-derivatives obtained via term-by-term differentiation. Since we supposedly can determine everything explicitly, we can expect this to be quite easy.
However it's not clear whether this actually converges to a smooth function, or if term-by-term differentiation is valid. Even in simple cases where $\phi_i,\psi_j$ are trigonometric and we get a Fourier series, we need very strong decay assumptions to get uniform convergence of $u$ and its derivatives (in general this need not be true).
This is where elliptic theory is useful, namely the estimates. Under appropriate assumptions, one can prove an estimate of the form,
$$ \lVert u \rVert_{C^2(\bar\Omega)} \leq C\left( \lVert f \rVert_{H^m(\Omega)} + \lVert u \rVert_{L^2(\Omega)}\right). $$
This is super useful, because now we know our solution is actually regular and we can differentiate it term-by-term. Of course this doesn't give us everything, since this doesn't prove the series converges to $u$ in a nice way; just that this $L^2$ limit is $C^2.$ If we want to prove further convergence/approximation results, some more work is necessary.