Let $\mathbb{B}$ be the unit ball in $\mathbb{R}^n$. Then it is true that if $u\in L^1(\mathbb{B})$ such that $\Delta(u)=0$ on $\mathbb{B}\backslash\{0\}$, then $\Delta(u)$, as a distribution on $\mathbb{B}$, is supported at the origin and has order at most one. This can be shown by noting that the fundamental solution to $\Delta(F)=\delta_{0}$ is a constant multiple of $\log|x|$ (for $n=2$) and $1/|x|^{n-2}$ (for $n\geq 3$) so its second derivatives do not belong to $L^1$.

Does this remain true for general elliptic differential operators? More specifically, let $D$ be an elliptic differential operators of order $N$ with smooth coefficients on $\mathbb{B}$. Suppose $u\in L^1(\mathbb{B})$ such that $D(u)=0$ on $\mathbb{B}\backslash\{0\}$. Consider $D(u)$ as a distribution on $\mathbb{B}$. It is obvious that the order of $D(u)$ is at most $N$ but must it always be strictly smaller than $N$? What I need is the case $D$ having polynomial coefficients with order $N=n+2$ so if the conclusion is not true in general, would it be true in this case?

Update (October 14, 2020): while I still cannot answer the question as originally stated, I think I have a proof if $D$ has smooth coefficients on $\mathbb{R}^n$ and $D(u)=0$ on all $\mathbb{R}^n\backslash\{0\}$. Then $D(u)=L(\delta)$ for some differential operator $L$ with constant coefficients and $\delta$ is the Dirac distribution at the origin. My idea is to use the Fourier transform. Indeed, the Fourier transform $\widehat{D(u)}(\zeta)$, since $u$ is compactly supported and $D$ has order $N$, is $o(|\zeta|^{N})$. To see this, note that $$\widehat{D(u)}(\zeta) = \int u(x)D^{*}(e^{-i\zeta\cdot x})dx = \sum_{|\alpha|\leq N}\zeta^{\alpha}\int u(x)c_{\alpha}(x)e^{-i\zeta\cdot x}dx$$ for some smooth $c_{\alpha}(x)$. These functions come from $D^{*}$. Since $u$ has compact support, Riemann-Lebesgue's Lemma shows that $\int u(x)c_{\alpha}(x)e^{-i\zeta\cdot x}dx\rightarrow 0$ as $|\zeta|\rightarrow 0$. So $\widehat{D(u)}(\zeta) = o(|\zeta|^{N})$. On the other hand, $\widehat{L(\delta)}(\zeta)$ cannot be $o(|\zeta|^{N})$ unless $L$ has order at most $N-1$. Therefore, $D(u)$ must have order at most $N-1$. The argument can be extended with some additional work to the case $D(u)=0$ on all $\mathbb{R}^n$ except at finitely many points. Note that this argument works without the assumption that $D$ is elliptic.

Updated (October 16, 2020): using the above argument and suggestion by Giorgio Metafune (see comments below), I believe we can show that the answer to my original question is affirmative. We don't need to assume that $D$ is elliptic as long as $u$ is smooth enough on $\mathbb{B}\backslash\{0\}$ so $D(u)(x)=0$ in the classical sense for $x\neq 0$. Since the proof seems rather standard (Fourier transform and cut-off functions), I would suspect that this should be known in the literature. If anyone has such a reference, please let me know.