…in R of course! There is a handy function to do those calculations. Normally (ahh!) you might resolve to a symbolic calculation package (Maple,Mathematica etc.) but that is not the situation any more. The calculations are done with the mnormt package. Relevant functions exist in other packages as well (R : Distributions) x <- seq(-2,4,length=21) y <- 2*x+10 z