Consider the following time-radial partial differential equation:
$$U_t(t,r) = \nabla^2 U(t,r) = U_{rr}(t,r) + {1\over r}U_r(t,r) \tag 1$$
If we would want to solve the above PDE numerically by using the finite difference method, the singularity at $r = 0$ would be a problem. For symmetric problems, the singularity point can be bypassed using the Maclaurin expansion of $U_r(t,r)$:
$$ U_r(t,r) \approx U_r(t,0)+r U_{rr}(t,0) \tag 2$$
Because the problem is symmetric, $U_r(t,0)$ is equal to zero. So instead of solving equation $(1)$ at $r=0$, we can solve:
$$U_t(t,0) = 2U_{rr}(t,0) \tag 3$$
However, I don't know how to generalize this approach to solving higher-order radial Laplacian equations. In my particular case, I want to solve the following PDE using the finite difference method:
$$U_t(t,r) = \nabla^4 U(t,r) + \nabla^2 U(t,r) = U_{rrrr}(t,r) + {2 \over r} U_{rrr}(t,r) - {1 \over r^2}U_{rr}(t,r) + {1 \over r^3}U_{r}(t,r) + U_{rr}(t,r) + {1\over r}U_{r}(t,r) \tag 4$$
Here is what I tried. I did the following Maclaurin expansions:
$$U_{rrr}(t,r)=U_{rrr}(t,0) + rU_{rrrr}(t,0)\tag 5$$ $$U_{rr}(t,r)=U_{rr}(t,0) + rU_{rrr}(t,0)+{r^2 \over 2} U_{rrrr}(t,0) \tag 6$$ $$U_r(t,r)=U_r(t,0) + rU_{rr}(t,0)+{r^2 \over 2} U_{rrr}(t,0) + {r^3 \over 6} U_{rrrr}(t,0) \tag 7$$
By inserting $(5)$, $(6)$ and $(7)$ into $\nabla^4U(t,r)$ I obtained:
$$ \nabla^4U(t,r) \approx \Bigg(2+ {1 \over 6}\Bigg)U_{rrrr}(t,0) + \Bigg(1+{1 \over 2}\Bigg){1\over r}U_{rrr}(t,0) \tag 8 $$
As you can see, I didn't get rid of $1/r$ which is a problem at $r = 0$. Also, I am not sure if it is still compatible to approximate $\nabla^2U(t,0)$ with expression $(2)$ or do I need to add now higher order terms of the expansion.
How can I bypass the singularity problem I encountered while trying to solve equation $(4)$ with the finite difference method?