Sorry for being late. The answer is simple. The cumulative distribution function of Matlab ... is not the same as the cumulative distribution function defined in Bensimhoun. Let us take a look at the documentation of Matlab. They say:
"The multivariate normal cumulative probability at X is defined as the probability that a random vector V, distributed as multivariate normal, will fall within the semi-infinite rectangle with upper limits defined by X.for example, Pr{V(1) ≤ X(1),V(2) ≤ X(2),...,V(d) ≤ X(d)}.
y = mvncdf(xl,xu,mu,SIGMA) returns the multivariate normal cumulative probability evaluated over the rectangle with lower and upper limits defined by xl and xu, respectively."
This is clearly different from the cumulative distribution function as defined in Bensimhoun. Here, you have the probability that a random vector V, distributed as multivariate normal, will fall within the Mahanalobis distance R from the mean, that is, inside the ELLIPSOID defined by the mean and covariance with scaling factor R.
For example, in 1D, the cumulative function of Matlab coincides with the usual cumulative function $F$ of the normal distribution (see definition in p. 1 of Bensimhoun), while the "mahalanobis cumulative" function in Bensimhoun coincides with $\int_{\mu-\sigma x}^{\mu+\sigma x} \varphi(t)dt$ (see Bensimhoun, beginning of p. 2, and how the two cumulative functions are related in the 1D-case, but not in higher dimension). The n-D dimensional cases are the exact generalizations of the 1-D dimensional case.