disp.3d <- function(model, data, model.type, steps = 30, lims = {}) {
      # Get a bounding interval for data[,3]
      # This will be taken as our third dimension
   if (length(lims) == 6) {
      z.min <- lims[5]
      z.max <- lims[6]
   } else {
      z.min <- floor(min(data[,3]))
      z.max <- ceiling(max(data[,3]))
   }
     # Break the interval into a number of subintervals
   z.int <- seq(z.min, z.max, ((z.max - z.min)/steps))
     # For the other two variables we create a grid
   if (length(lims) == 6) {
      res <- f.create.grid(data.frame(data[,1], data[,2]), c(steps, steps), lims[1:4])
   } else {
      res <- f.create.grid(data.frame(data[,1], data[,2]), c(steps, steps))
   }

   xp <- res$xp
   yp <- res$yp
   gp <- res$zp      # The grid points
       # In order to display the separators we look for
       # changes in predicted class
   x.change <- {}
   y.change <- {}
   z.change <- {}
   Z.all <- {}
   for (i in 1:steps){
         # We will look at a grids on a sequence of planes
         # in the d.flea[,4] direction.
      data.3d <- cbind(gp, z.int[i])
         # The predict method likes the names to be the same
      colnames(data.3d) <- colnames(data)
         # We use the model to predict the class of every point
         # on the grid on the current data[,3] plane
      if (model.type == "lda") {
         Z <- unclass(predict(model, data.3d)$class)
       } else if (model.type == "qda") {
         Z <- unclass(predict(model, data.3d)$class)
       } else if (model.type == "net") {
         Z <- max.col(predict(model, data.3d))
       }
       Z.all <- rbind(Z.all, cbind(data.3d, as.numeric(Z)))
         # Put the result into a matrix
      Z.m <- matrix(Z, nrow = steps + 1)
         # Find the positions in a row where the class changes
      Z.m.diff.r <- t(apply(Z.m, 1, diff))
      cc <- which(Z.m.diff.r != 0, arr.ind = T)
         # Find the positions in a column where the class changes
      Z.m.diff.c <- apply(Z.m, 2, diff)
      rr <- which(Z.m.diff.c != 0, arr.ind = T)
         # Convert the matrix positions to grid positions
      x.change <- c(x.change, xp[rr[,1]],xp[cc[,1]])
      y.change <- c(y.change, yp[rr[,2]],yp[cc[,2]])
      z.change <- c(z.change, rep(z.int[i], length(rr[,1]) + length(cc[,1])))
   }
      # After we have found the change points on all the grids
   list(changePts = data.frame(x.change, y.change, z.change), classes = Z.all)
}