R help - Nov 2009 - Partial derivatives of the multivariate cumulative distribution

I'm currently using the mvtnorm package to model unobserved  
heterogeneity in a structural model and using optim to estimate the  
model. I have got good clues that convergence is not really a problem  
but the hessian matrix estimate is very bad. To overcome this problem,  
I'm constructing an OPG estimator of the information matrix and I was  
wondering if there were an easy way to obtain partial derivatives of  
say for instance:

P1  <-  
pmvnorm(lower=c(-Inf,-Inf,-Inf,-Inf),upper=c(theta1,theta2,theta3,theta4),corr=ssigma)

with respect to the mean parameters theta1, theta2, theta3, theta4 and  
the non-diagonal parameters in sigma, hence $\partial P_1 / \partial  
\theta_1$, etc...

I can deal with numerical or analytical partial derivatives - a  
gradient would be fine since all observations share the same partial  
derivative.

Stephane