my.bp <- function(formula,
                  data) {

  data.f <-  model.frame(formula = formula,
                         data = data)

  # Get all the basic IV and DV info
  
  iv.name <- names(data.f)[2]
  dv.name <- names(data.f)[1]

  iv.vec <- data.f[,2]
  dv.vec <- data.f[,1]

  iv.levels <- sort(unique(iv.vec))

  n.iv <- length(iv.levels)

  # Set up plotting space
  
  plot(1,
       xlim = c(0, n.iv + 1),
       ylim = c(min(dv.vec), max(dv.vec)),
       type = "n",
       ylab = dv.name,
       xlab = iv.name,
       xaxt = "n"
  )

  # Loop over all iv.levels
  
  for(i in 1:n.iv) {

    # Get dv for current iv level

    dv.current <- dv.vec[iv.vec == iv.levels[i]]

    # Add plotting element for current iv level

    boxplot(dv.current,
            at = i,
            add = T
    )



  }

}

Testing

my.bp(formula = weight ~ Time,
      data = ChickWeight)