Dear Rers, It seems the usual sum function can work. Anyway, I appreciate your time on this. Best wishes, Liqiu Dear Rers, I am trying to do a 2 dimmensional intergration for a function. the function is summation of another function evaluated at a series of vector values. I have difficulty to code this function. For example: I have function f which is a bivariate normal density function: #define some constants err<-0.5 m<-5 times<-seq(0, m-1) rou<-sum(times)/sqrt(m*sum(times^2)) sig.w<- sqrt(m*err) sig.wt<-sqrt(sum(times^2)*err) #bivariate normal density f<-function(x, y, u.x, u.y) exp(-((x-u.x)^2/sig.w^2+(y-u.y)^2/sig.wt^2-2*rou*(x-u.x)*(y-u.y)/(sig.w*sig.wt))/(2*(1-rou^2)))/(2*pi*sig.w*sig.wt*sqrt(1-rou^2)) ### I would like to have a function g which is defined as ###### uw = 1:n uwt = (n+1):2n g = function(x, y) f(x,y, uw[1], uw[1])+f(x,y, uw[2], uwt[2])+ ...+f(x,y, uw[n], uwt[n]) ####### if n is very large, I am not able to write all them down, How can I code the function g. Thank you for your consideration. Best wishes, Liqiu