Consider Schrodinger equation (where $E, \mu$ and $\hbar$ are constants and the form of $V(r)$ was given)
$$E \psi (\vec x) = \frac{-\hbar^2}{2\mu} \nabla^2 \psi (\vec x) + V(r) \psi (\vec x)$$
or in a more expanded form in spherical coordinates
$$E\psi(\vec x) = -\frac{\hbar^2}{2\mu r^2}\frac{\partial}{\partial r}\bigg(r^2\frac{\partial\psi(\vec x)}{\partial r}\bigg) + \frac{1}{2\mu r^2}\vec L ^2\psi(\vec x ) + V(r)\psi(\vec x) $$
where we have
$$\vec L^2 =-\hbar^2\bigg[\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\bigg(\sin\theta\frac{\partial}{\partial\theta}\bigg) + \frac{1}{\sin^2\theta} \frac{\partial^2}{\partial^2\phi}\bigg]$$
I read a book that said
if $V(r)$ is not extremely singular at $r=0$ then wave function $\psi$ must be a smooth function of the Cartesian components $x_i$ near $x = 0$, in the sense that it can be expressed as a power series in these components.
So I have two questions about this statement
1. What is mathematical theorem or something else, ensures that we can do this and expand the $\psi(\vec x)$ in power series in near of center?
2. What is the boundary between extremely singular and not extremely singular in this situation?
Partial Answer
Assuming you use separation of variables to solve for $\psi(\vec{r})=R(r)\cdot Y^m_l(\theta,\phi)$, where the $Y^m_l(\theta,\phi)$ are spherical harmonics, which are eigenfunctions of $L^2$, then I believe the Frobenius method is useful. (I've seen the method presented as a theorem in ODE books like Simmons).
Essentially when solving an ODE of the form:
$$y''+P(x)y'+Q(x)y=0$$
where one or both of $P(x)$ and $Q(x)$ are singular (non-analytic) at a point $x=0$ say, with $xP(x)$ and $x^2Q(x)$ analytic at $x=0$, then there is a solution of the form
$$x^m \sum_{n=0}^\infty a_n x^n$$
where $m$ is a root of $m(m-1)+mp_0+q_0=0$, with $p_0=\lim_{x\to0}xP(x)$ and $q_0=\lim_{x\to0}x^2Q(x)$. The second solution depends on whether there are repeated roots to this equation or whether the difference between the roots is an integer.
To bring this back to the Schrodinger equation, with $L^2 Y^m_l(\theta,\phi)=\hbar^2 l(l+1)$ you can separate and solve the radial equation
$$-\lambda R(r)=\frac{d^2 R(r)}{dr^2}+\frac{2}{r}\frac{d R(r)}{dr}-\left(\frac{l(l+1)}{r^2}+\frac{2\mu V(r)}{\hbar^2}\right) R(r)$$
with $\lambda=\frac{2\mu E}{\hbar^2 }$. This is of the above ODE form for the frobenius method with $P(r)=\frac{2}{r}$ and $Q(r)=\lambda-\frac{l(l+1)}{r^2}-\frac{2\mu V(r)}{\hbar^2}$.
So I believe to use the frobenius method you want $V(r)$ diverging no worse than $\frac{1}{r^2}$ at $r=0$, although you might need extra conditions to guarantee $\psi$ is defined at $r=0$, or to have some restrictions on the allowable energies (minimum energy for ground state - energies bounded below).
Update
If you want to see that you can separate the variables and don't know of the spherical harmonics (yet), you can see that the $L^2$ operator only acts on the angular parts, and this is being added to terms that only operate on the radial part. So you can ask for solutions of the form $\psi(\vec{r})=R(r)\cdot \chi(\theta,\phi)$, subbing this into the Schrodinger equation you get
$$ER(r)\chi(\theta,\phi)= -\frac{\hbar^2}{2\mu}\left(\frac{d^2}{d r^2}R(r)\right)\chi(\theta,\phi) -\frac{\hbar^2}{\mu r}\left(\frac{d}{d r}R(r)\right)\chi(\theta,\phi)+ \frac{1}{2\mu r^2}\vec L ^2\left(\chi(\theta,\phi)\right) R(r) + V(r)R(r)\chi(\theta,\phi)$$
which you can rearrange and divide across by $R$ and $\chi$ (at points/regions where they're not zero) and multiply by $r^2$ to get
$$-\frac{1}{2\mu \chi(\theta,\phi)}\vec L ^2\left(\chi(\theta,\phi)\right)= -\frac{\hbar^2 r^2}{2\mu R(r)}\left(\frac{d^2}{d r^2}R(r)\right) -\frac{\hbar^2 r}{\mu R(r)}\left(\frac{d}{d r}R(r)\right)+ r^2(V(r)-E)$$
Now the right hand side is independent of $\theta$ and $\phi$ and the left hand side is independent of $r$, so if you take partial derivates, you can conclude that both sides must be a constant independent of all of $\theta,\phi,r$, call it $A$. This gives you two separate equations to solve.
$$ \vec L ^2\left(\chi(\theta,\phi)\right)=-2\mu A \chi(\theta,\phi)$$ and $$-\frac{\hbar^2 r^2}{2\mu}\left(\frac{d^2}{d r^2}R(r)\right) -\frac{\hbar^2 r}{\mu}\left(\frac{d}{d r}R(r)\right)+ r^2(V(r)-E) R(r)=A R(r)$$