I'm trying to show $$\nabla \times \vec{L} = \frac{1}{i}(\vec{r}\nabla^{2}-\nabla(1+r\frac{\partial}{\partial{r}})) $$ where $ \vec{L}\psi = \frac{1}{i}(\vec{r}\times \nabla)\psi$.
I'm able to expand and simplify all the terms except for the $(\nabla\psi \cdot\nabla)\vec{r}$ term. From Jackson E&M, it seems that the term should equal to $\nabla\psi$ since it will give me the correct result.
In this answer, I denote the cartesian unit vectors by $\mathbf{e}_i$ ($i=1,2,3$) and the radius vector $\mathbf{r}=\sum_ix_i\mathbf{e}_i$.
Assuming $\mathbf{L}=\frac{1}{\mathrm{i}}\mathbf{r}\times\nabla$ we have \begin{align} \nabla\times\mathbf{L}&=\sum_{i,j,k}\varepsilon_{ijk}\mathbf{e}_i\partial_jL_k\\ &=\sum_{i,j,k}\varepsilon_{ijk}\mathbf{e}_i\partial_j\left(\frac{1}{\mathrm{i}}\mathbf{r}\times\nabla\right)_k\\ &=\frac{1}{\mathrm{i}}\sum_{i,j,k}\varepsilon_{ijk}\mathbf{e}_i\partial_j\sum_{l,m}\varepsilon_{klm}x_l\partial_m\\ &=\frac{1}{\mathrm{i}}\sum_{i,j,k,l,m}\varepsilon_{kij}\varepsilon_{klm}\mathbf{e}_i\partial_j\left(x_l\partial_m\right). \end{align} Where I used $\varepsilon_{ijk}=\varepsilon_{kij}$. Now \begin{equation} \sum_k\varepsilon_{kij}\varepsilon_{klm}=\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl} \end{equation} so \begin{align} \nabla\times\mathbf{L}&=\frac{1}{\mathrm{i}}\sum_{i,j,l,m}\left(\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl}\right)\mathbf{e}_i\partial_j\left(x_l\partial_m\right)\\ &=\frac{1}{\mathrm{i}}\sum_{i,j,l,m}\delta_{il}\delta_{jm}\mathbf{e}_i\partial_j\left(x_l\partial_m\right)-\frac{1}{\mathrm{i}}\sum_{i,j,l,m}\delta_{im}\delta_{jl}\mathbf{e}_i\partial_j\left(x_l\partial_m\right)\\ &=\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_i\partial_j\left(x_i\partial_j\right)-\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_i\partial_j\left(x_j\partial_i\right). \end{align} To rewrite the summands, let's evaluate their action on a twice continuously differentiable scalar field $f(x_1,x_2,x_3)$. We get \begin{align} \partial_j\left(x_j\partial_if\right)&=(\partial_jx_j)(\partial_if)+x_j\partial^2_{ij}f\\ &=\delta_{jj}\partial_if+x_j\partial^2_{ij}f. \end{align} Similarly \begin{align} \partial_j\left(x_i\partial_jf\right)&=(\partial_jx_i)(\partial_jf)+x_i\partial^2_{j}f\\ &=\delta_{ij}\partial_jf+x_i\partial^2_{j}f. \end{align} Since $f$ was arbitrary enough we have the operator equalities \begin{equation} \partial_j\left(x_i\partial_j\right)=\delta_{ij}\partial_j+x_i\partial^2_{j},\quad \partial_j\left(x_j\partial_i\right)=\delta_{jj}\partial_i+x_j\partial^2_{ij}. \end{equation} Hence \begin{align} \nabla\times\mathbf{L}&=\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_i\left(\delta_{ij}\partial_j+x_i\partial^2_{j}\right)-\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_i\left(\delta_{jj}\partial_i+x_j\partial^2_{ij}\right)\\ &=\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_i\delta_{ij}\partial_j+\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_ix_i\partial^2_{j}-\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_i\delta_{jj}\partial_i-\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_ix_j\partial^2_{ij}\\ &=\frac{1}{\mathrm{i}}\nabla+\frac{1}{\mathrm{i}}\mathbf{r}\nabla^2-\frac{3}{\mathrm{i}}\nabla-\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_ix_j\partial^2_{ij}\\ &=-\frac{2}{\mathrm{i}}\nabla+\frac{1}{\mathrm{i}}\mathbf{r}\nabla^2-\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_i\partial_i(x_j\partial_j)+\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_i\underbrace{\partial_i(x_j)}_{=\delta_{ij}}\partial_j\\ &=-\frac{2}{\mathrm{i}}\nabla+\frac{1}{\mathrm{i}}\mathbf{r}\nabla^2-\frac{1}{\mathrm{i}}\sum_{i,j}\mathbf{e}_i\partial_i(x_j\partial_j)+\frac{1}{\mathrm{i}}\sum_{i}\mathbf{e}_i\partial_i\\ &=-\frac{2}{\mathrm{i}}\nabla+\frac{1}{\mathrm{i}}\mathbf{r}\nabla^2-\frac{1}{\mathrm{i}}\nabla(\mathbf{r}\cdot\nabla)+\frac{1}{\mathrm{i}}\nabla\\ &=\frac{1}{\mathrm{i}}\mathbf{r}\nabla^2-\frac{1}{\mathrm{i}}\nabla-\frac{1}{\mathrm{i}}\nabla(r\partial_r)\\ &=\frac{1}{\mathrm{i}}\left[\mathbf{r}\nabla^2-\nabla(1+r\partial_r)\right],\quad\text{Q.E.D.} \end{align}
Indeed, your identity $$ (\nabla \psi)\cdot (\nabla \vec r) = (\nabla \psi ) $$ holds, where this parenthesis indicates the gradient only acts on the argument (ψ) in it, a dangerously ambiguous notation!
Further note your expression $$ \nabla ~~ r\partial_r = \nabla ~~\vec r \cdot \nabla = (1+\vec r\cdot \nabla ) \nabla ~. $$ So, your rectified expression trivially coincides with the simpler result $$ \vec A= \nabla \times (\vec r \times \nabla) =\vec r \nabla^2 -(2+ \vec r\cdot \nabla) \nabla ~ \tag 1 $$ where, as an easy check, you observe that acting on the function $r^2$ you get zero, $\vec A r^2=0$, as $r^2$ is in the kernel of the operator $\vec L$.
To prove (1), you just note $$ \nabla \cdot \vec r= \vec r\cdot \nabla + 3, $$ so that $$ \nabla \times (\vec r \times \nabla) = \nabla + \vec r \nabla^2 - \nabla \cdot \vec r ~~ \nabla= \nabla + \vec r \nabla^2 -\vec r \cdot \nabla -3\nabla . $$
Your operator $\vec A$ is "half" the LRL vector.
Writing $\nabla \vec r$ is pretty horrible (as others comment) but if you choose a better notation and write instead $\partial_a x^b = \delta^b_a$ I think you will find it is not difficult.
