I am calculating the derivative of the following function with respect to $\boldsymbol{x}$, where $\boldsymbol{x}$ is a vector of length n and $\phi$ denotes the multivariate standard normal distribution:
$$f=\int_{-\infty}^{t_1=x_1}\int_{-\infty}^{t_2=x_2}...\int_{-\infty}^{t_n=x_n} \boldsymbol{t} \phi(\boldsymbol{t}) d\boldsymbol{t}$$
I thought the answer was this
$$\frac{d f}{d\boldsymbol{x}}=\sum_{i=1}^n g_i(x_i),$$
where $$ g_i(x_i)=\int_{-\infty}^{t_1=x_1}...\int_{-\infty}^{t_{i-1}=x_{i-1}}\int_{-\infty}^{t_{i+1}=x_{i+1}}...\int_{-\infty}^{t_n=x_n} [t_1... t_{i-1}x_it_{i+1}...t_n] \phi([t_1... t_{i-1}x_it_{i+1}...t_n]) dt$$
However, when I evaluated this result numerically in R, it was incorrect. Is my solution wrong?