There’s a relatively new R package that has code to calculate he cdf of a multivariate normal that uses the powell dogleg method TruncatedNormal/mvNcdf.R at master · lbelzile/TruncatedNormal · GitHub. Apparently if that fails it tries a convex solver. Anyway, the method calculates the jacobian (TruncatedNormal/gradpsi.R at 6a5428e3861e8ff77ac55cf9c1526dfcb3d9ce4a · lbelzile/TruncatedNormal · GitHub) and passes that to the solver.

@charlesm93 is it possible to give the solver the jacobian in Stan if I wanted to code this up? Or do I have to go to c++?