I am basically checking that the answer which Mathematica gives for computing twice the grad of a function is the same as the one which I compute explicitly myself. So I define the function mu_1 and apply grad twice:
\[Phi] = 1/(4*pi*R);
\[Mu]1 = \[Phi]*Exp[-gam1*R];
\[Mu]11 = FullSimplify[Grad[\[Mu]1, r]];
\[Mu]12 = FullSimplify[Grad[\[Mu]11, r]];
However, when I try the explicit computation by applying the formula for grad in spherical coordinates I get the following:
$\mu_{\gamma}=-\frac{e^{-\gamma}}{4 \pi r}$
$\nabla \mu_{\gamma} = \frac{\textbf{r}}{|\textbf{r}|} \frac{\text{d}\mu_{\gamma}}{\text{d}|\textbf{r}|}$.
$\nabla (\nabla \mu_{\gamma}) = \frac{\nabla (\textbf{r})}{|\textbf{r}|}\frac{\partial}{\partial r}\mu_{\gamma} + \hat{r} \nabla \bigg( \frac{\partial}{\partial r} \mu_{\gamma} \bigg)=\frac{\textbf{r}}{|\textbf{r}|} \frac{\text{d} \mu_{\gamma}}{\text{d} |\textbf{r}|} + \frac{\textbf{r}}{|\textbf{r}|}\frac{\textbf{r}}{|\textbf{r}|} \frac{\partial }{\partial r} \bigg( \frac{\partial}{\partial r} \mu_{\gamma} \bigg) = \frac{\delta_{ij}}{|\textbf{r}|}\frac{\text{d} \mu_{\gamma}}{\text{d}|\textbf{r}|} + \frac{\textbf{r} \otimes \textbf{r}}{|\textbf{r}|^2} \frac{\text{d}^2 \mu_{\gamma}}{\text{d} |\textbf{r}|^2}. $
After performing the derivatives to get the expression, if I then try to subtract the result from Mathematica from the result of my brief computation, I do not get 0, so cl