How can I add contour lines to a 3d plot built with persp()?

I am trying to plot a 3D chart of the function x^2+y^2 and add the contour lines to try to have the same aspect as this figure built with Matlab:

However, I am not able to add the contour lines, the grid lines and the math notation in the axis. This is the plot I got with the code provided using the persp function:

library(rgl)
library(rsm)

# Generate data for the plot
x <- seq(-5, 5, length.out = 25)
y <- seq(-5, 5, length.out = 25)
z <- outer(x, y, function(x, y) x^2 + y^2)

jet.colors <- colorRampPalette( c("blue", "red") ) 
pal <- jet.colors(100)
col.ind <- cut(z,100) 

# Create the plot
persp(x, y, z, 
zlab = "Q1(x1,x2)",
xlab = expression(x^2), #not working
ylab = "x1",
theta = 50, phi = 20, expand = 0.5, 
ticktype = "detailed",
contour(x, y, z),
col = pal[col.ind])

Your 2nd last argument to persp, i.e. contour(x, y, z), is being interpreted as an xlim argument, because you didn't give it a name. As far as I know persp() doesn't support adding 2D annotations in the same call.

Your code specifies library(rgl) and library(rsm) but as far as I can see you don't use either of those packages.

To do what you want to do, I can see two approaches. The first approach is to draw the plot with the correction mentioned above, then use the trans3d function to draw contour lines on it. Since there's no hidden line removal, you'd need to draw the main plot again to cover the lines that should have been hidden.

Regarding the question about your axis label x^2: see https://stackoverflow.com/a/23420409/2554330 for a solution using lattice graphics, or figure out the coordinates you want, and use text(trans3d(x,y,z, trans), label=expression(x^2)) to put it there.

Here's an attempt at all of this:

# Generate data for the plot
x <- seq(-5, 5, length.out = 25)
y <- seq(-5, 5, length.out = 25)
z <- outer(x, y, function(x, y) x^2 + y^2)

jet.colors <- colorRampPalette( c("blue", "red") ) 
pal <- jet.colors(100)
col.ind <- cut(z,100) 

# Create the plot and save the transformation
trans <- persp(x, y, z, 
      zlab = "Q1(x1,x2)",
      xlab = "",
      ylab = "x1",
      theta = 50, phi = 20, expand = 0.5, 
      ticktype = "detailed",
      col = pal[col.ind])

# Compute the contour lines and add them
lines <- contourLines(x, y, z)
bottom <- min(z)
for (i in seq_along(lines)) {
  segment <- lines[[i]]
  lines(trans3d(segment$x, segment$y, bottom, trans))
}

# The lines overwrite the surface, so redraw it
par(new = TRUE)  # to allow us to draw to the same plot
persp(x, y, z, 
      zlab = "Q1(x1,x2)",
      xlab = "", #not working
      ylab = "x1",
      theta = 50, phi = 20, expand = 0.5, 
      ticktype = "detailed",
      col = pal[col.ind])

# Add the xlab
text(trans3d(0, -6, bottom - 10, trans), labels=expression(x^2))

^{Created on 2024-01-29 with reprex v2.0.2}

You could also add grid lines if you want before redrawing the plot (or after, if you want them in front of it).

Another approach uses rgl. It has a function show2d() that draws a 2D plot on a rectangle in a 3D scene.