Axis labeling for subplots
Setting custom axis labels for 3d plots doesn't work with Plots and plotlyjs() backend. Only with GR backend your labels are displayed.
You can try this version using PLotlyJS.jl instead Plots.jl:
fig=make_subplots(rows=1, cols=2, specs =[Spec(kind="scene") Spec(kind="xy")], 
                  horizontal_spacing=-0.1, column_widths=[0.65, 0.35])
add_trace!(fig, PlotlyJS.surface(x=sol.t, y=collect(x), z=z', showscale=false), row=1, col=1)
add_trace!(fig, PlotlyJS.contour(x=sol.t, y=collect(x), z=z), row=1, col=2)
relayout!(fig, template=templates["plotly_white"], font_size=11,
                width=1000, height=600, scene=attr(xaxis_title="Time", yaxis_title="Space", 
                zaxis_title="|u|²", camera_eye=attr(x=1.8, y=1.8, z=1)),
                xaxis2_title="Time", yaxis2_title="Space",  margin_l=15) 
display(fig)



The following code generates this image 
I want the "y"-axis label to be "Space" and the "x"-axis label to be "time" for the left subplot. However, I am failing to achieve this. Why does my plotting code not do as I desire?
p1 = surface(sol.t, x, z, xlabel="Time", ylabel="Space", zlabel="|u|²", colorbar = false)
p2 = contourf(sol.t,x,z, xlabel="Time", ylabel="Space")
plt = plot(p1,p2,layout=(1,2), size=(1200,800))


using DifferentialEquations, LinearAlgebra, Plots, SparseArrays
plotlyjs()

N₁=31 # Number of waveguides / size of solution vector
γ=1  # Nonlinear term strength parameter
h=1 # Grid spacing 

centerGrid  = (N₁-1)/2;
x = -centerGrid:centerGrid;

# Coefficient matrix of second-order centered-difference operator (δ²u)ₙ
M           = spdiagm(-1 => fill(1,N₁-1), 0 => fill(-2,N₁), 1 => fill(1,N₁-1))
M[N₁,1]     = 1; # Periodic boundary conditions
M[1,N₁]     = 1;

# RHS of DNLS. The solution vector u is a N₁x1 complex vector
g₁(u,p,t)   = 1*im*(p[1]*M*u + @.(γ*((abs(u))^2).*u) )

# Julia is explicitly typed (e.g, cannot have Int and Complex in same array) and so we must convert the object containing the initial data to be complex
u0  = Complex.(sech.(x))

tspan = (0.0,200)
prob = ODEProblem(g₁,u0,tspan, [h])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

z= [abs(sol.u[i][j])^2  for j=1:N₁, i=1:size(sol)[2]] # |u|²

p1 = surface(sol.t, x, z, xlabel="Time", ylabel="Space", zlabel="|u|²", colorbar = false)
p2 = contourf(sol.t,x,z, xlabel="Time", ylabel="Space")
plt = plot(p1,p2,layout=(1,2), size=(1200,800))

