Thank you very much for your insightful comment. Sometimes EJMR works after all.

As to your comment about approximations for any model that I worked until now analytical derivatives were considerably faster yet as you mentioned was hard to derive. Of course, it might be the case that until now I haven't estimated things as nasty as you might have done. But I can say that I estimated non-traditional MLs with non-normal error terms.

So, why don't you estimate a random effects model instead? I could not get what you mean by interval restrictions though. Be careful about the standard errors that you are getting. I am not sure whether the inverse Hessian works when there are restrictions in the model.

Random effects models are even harder to estimate. Interval restrictions are just constraints on the parameter space. And yes: the standard errors if a constraint are binding is a total mess. Of course, you could use a soft constraint, but then your standard errors would mostly be telling you about the soft constraint.

If you must use MLE then I would strongly recommend you finding gradient analytically which would boost the speed and accuracy of the estimates. If possible find the Hessian as well. Again this will boost the speed and increase the accuracy further. Good luck.

Read up on CG as mentioned above. So long as you can find an approximation to the gradient which converges locally, you do not need it analytically or even exactly. Inverse Hessian can be found through approximations (SR1 or DFP2). Sometimes, if the approximations are close and the analytical forms are nasty, the approximation methods can even be faster.