Hello David and Dennis, thank you two for your replies. You not only helped me solve my problem, but also helper me to understand my problem better. Here is my -now working- code: x<-y<-seq(-4,4,len=40) g<-function(a,b) { dmvnorm(x=cbind(a,b),sigma=matrix(c(4,2,2,3), ncol = 2)) } z<-outer(x,y,g) persp(x,y,z) Thank you a lot! Severin