Integration by parts to prove a multivariable identity Let $X$ denote a normal random variable and $f$ some differentiable function.
In single-variable calculus, I can prove that :
$\mathbb{E}[Xf(X)]=\mathbb{E}[X^2]\mathbb{E}[f'(X)]$.
by integrating by parts.
Indeed, $\mathbb{E}[Xf(X)]=\int_{-\infty}^{+\infty}xf(x)w(x)dx$ with $w(x)=\frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac12 \frac{x^2}{\sigma^2}}$
so it suffices to notice that $\frac{dw(x)}{dx}=-\frac{x}{\sigma^2}w(x)$ to complete the integration by parts and obtain :
$\mathbb{E}[Xf(X)]=\sigma^2\int_{-\infty}^{+\infty}f'(x)w(x)dx=\mathbb{E}[X^2]\mathbb{E}[f'(X)]$
Now, some exercise says it can be generalized to vectors :
$\mathbb{E}[X_if(X_1,...,X_n)]=\Sigma_{j=1}^{n}\mathbb{E}[X_iX_j]\mathbb{E}[\frac{\partial f}{\partial X_j}(X_1,...X_n)]$ $\space \space \space$(*)
As far as I know, the probability density function for multinormal random variables is usually defined as $w(\vec X)=Ne^{-\vec x^TC\vec x}$ with $N$ the normalization factor and C the covariance matrix whose inverse yields the $\sigma_{ij}$.
It's been far too long since the last time I studied multivariable calculus (which was barely formally taught to me actually, $2$-$3$ maths  exercises max for my whole life as a graduate physics student, rest was just knowing by heart $r^2sin(\theta)$ without ever needing the formal definition of a Jacobian), so now I am having trouble...
What I have tried :
I first notice  that in the RHS of (*), the vector differential $\frac{\partial}{\partial \vec Y}$ is involved on $f(\vec X)$, with $\vec Y=(c_{i1}X_1,c_{i2}X_2,...,x_{in}X_n)$ (should be a column vector actually)
Naively generalizing integration by parts, I get the following :
$\int_{-\infty}^{+\infty}\frac{\partial f(\vec X)}{\partial \vec Y}w(\vec X)d\vec X=\color{blue}{\int_{-\infty}^{+\infty}\frac{\partial (f(\vec X)w(\vec X))}{\partial \vec Y}d\vec X} - \int_{-\infty}^{+\infty}f(\vec X)\frac{\partial w(\vec X)}{\partial \vec Y}d\vec X$
If I remember correctly, $\frac{d\vec X}{d \vec Y}$ is the Jacobian matrix, and since $\vec Y$ is just $\vec X$ times a diagonal matrix with diagonal entries $(c_{ij})_{j=1,...,n}$, I have that $\frac{d\vec X}{d \vec Y}=(c_{ij}^{-1})_{j=1,...,n}=(\sigma_{ij})_{j=1,...,n}$
So the blue part is just $\sum_{i=1}^{n} \sigma_{ij}[f(X_i)w(X_i)]_{-\infty}^{+\infty}=0$
Therefore, I am left with :
$\int_{-\infty}^{+\infty}\frac{\partial f(\vec X)}{\partial \vec Y}w(\vec X)d\vec X=- \int_{-\infty}^{+\infty}f(\vec X)\frac{\partial w(\vec X)}{\partial \vec Y}d\vec X$
and I have to prove that : 
$- \int_{-\infty}^{+\infty}f(\vec X)\frac{\partial w(\vec X)}{\partial \vec Y}d\vec X = \mathbb{E}[X_if(X_1,...,X_n)] = \int_{-\infty}^{+\infty}X_if(\vec X)w(\vec X)d \vec X$
which could be done if : 
$-\frac{\partial w(\vec X)}{\partial \vec Y} = X_iw(\vec X)$
Since $w(\vec X)$ is supposed to be a scalar, I get that :
$-\frac{\partial w(\vec X)}{\partial \vec Y} =-\sum_{j=1}^n \sigma_{ij} \frac{dw(\vec X)}{dX_j}$
Using $\frac{dw(\vec X)}{dX_j}=-w(\vec X)\sum_{i=1}^n c_{ij}X_i$
I obtain :
$-\frac{\partial w(\vec X)}{\partial \vec Y} =w(\vec X)\sum_{k=1}^{n}\sum_{j=1}^n \frac{\sigma_{ij}}{\sigma_{kj}}X_k=w(\vec X)\sum_{k=1}^n X_k \sum_{j=1}^{n}\frac{\sigma_{ij}}{\sigma_{kj}}$
$=w(\vec X)X_i+\color{green}{w(\vec X) \sum_{k=1, k\neq i}^n X_k \sum_{j=1}^{n}\frac{\sigma_{ij}}{\sigma_{kj}}}$
which isn't what I am looking for since the green part doesn't seem to be $0$... and doesn't seem to vanish either when put into the integral as we don't know anything about $f$.
Can someone pinpoint where I am not using multivariable calculus correctly ?
 A: Sorry I am super tired and not completely following your notation so I can't exactly tell you just now where what you're doing is going wrong, but if I offer some pointers as to how this result could be derived which I hope will help. 
Note that if $\varphi(\textbf{x})$ is the density of a multivariate normal with covariance matrix $\Sigma$ then from the logarithmic derivative we have
$$
\nabla \varphi(\textbf{x}) = - \left(\Sigma^{-1}\textbf{x}\right)\varphi(\textbf{x}),
$$
now let $\Lambda = \Sigma^{-1}$, so that an individual component is just given by
$$
\frac{\partial \varphi}{\partial x_i} = - \left(\sum_j \Lambda_{ij} x_j \right) \varphi(\textbf{x}),
$$
now taking arbitrary $f$ and assuming appropriate boundary conditions we get
$$
\begin{align}
\int_{\mathbb{R}^n} \frac{\partial f(\textbf{x}) }{\partial x_i} \varphi(\textbf{x})dx &=\int_{\mathbb{R}^{n-1}} \left[\int\frac{\partial f}{\partial x_i} \varphi(\textbf{x}) dx_i \right]d\textbf{x}_{-i} \\
&=\int_{\mathbb{R}^{n-1}}\left[ - \int f(\textbf{x})\frac{\partial \varphi(\textbf{x})}{\partial x_i}\right] d\textbf{x}_{-i} \\
&= \int_{\mathbb{R}^{n-1}}\left[ \int f(\textbf{x})\sum_{j} \Lambda_{ij}x_j\varphi(\textbf{x})dx_i\right] d\textbf{x}_{-i} \\
&= \int_{\mathbb{R}^n} f(\textbf{x}) \left(\sum_j \Lambda_{ij} x_j \right) \varphi(\textbf{x})  d\textbf{x}
\end{align}
$$
So in terms of expectations and the inverse covariance (precision matrix) we get,
$$
\mathbb{E}\left[ \frac{\partial f}{\partial x_i} (\textbf{X}) \right]=\sum_{j}\Lambda_{ij} \mathbb{E}\left[ X_j \cdot f(\textbf{X})\right].
$$
