Take the following example.
library(dplyr)
temp <- data.frame(lapply(1:3, function(i) rnorm(5, 0, 1)))
names(temp) <- paste0("X", 1:3)
temp_each <-
    temp %>%
    mutate_each(funs(mean, median))
Examining the names of temp_each, we see that
> names(temp_each)
[1] "X1"        "X2"        "X3"        "X1_mean"   "X2_mean"   "X3_mean"   "X1_median" "X2_median" "X3_median"
that is, the final columns are in groups of three, always ordered  X1, X2, X3 + the function applied.
However, I would like it to look like this
[1] "X1"        "X1_mean"   "X1_median" "X2"        "X2_mean"   "X2_median" "X3"        "X3_mean"   "X3_median"
Does anyone know how to implement this, preferably using dplyr, for a data frame with many many columns and arbitrary column names?