Hi, is there a way of simplifying the following code: G <- rep(NA,n) for(i in 1:n) { gj <- 0 for(j in 1:n) { for(l in 1:n) { for(m in 1:n) { gj <- gj+G.fun(XB[i]+p[3]*X[j,3]+p[4]*X[l,4]+p[5]*X[m,5],ff) } } } G[i] <- gj/n^3 } Thanks. Joaquim Santos