How to create a Surface plot in Plot.jl from ```solve()``` data

MWE:

N₁=5
γ=-1 #nonlinear term strength parameter

g₁(u,p,t) = 1*im*(p[1]*Tridiagonal(fill(1,N₁), fill(-2,N₁+1), fill(1,N₁))*u + γ*u.^3)
vals=0.3*ones(N₁+1)
u0  =map(Complex,vals)


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

I want to create a surface plot of the absolute value squared data in sol similar to what is show here:

Attempt: I want to use the Plots.jl package and utilize the surface() function. The time data in sol can be accessed as sol.t and the spatial data can be accessed as sol.u. For example, at sol.t[1] the corresponding spatial data is a 6-component vector sol.u[2] given by 6-element Vector{ComplexF64}: 0.2999679677368571 - 0.0030866750511122525im 0.30001119699789686 - 0.000402829297042341im 0.2999991898889756 - 0.00040261260929592426im 0.2999991898889756 - 0.00040261260929592426im 0.30001119699789686 - 0.000402829297042341im 0.2999679677368571 - 0.0030866750511122525im

The components of the vector are the spatial data at the first, second, third, fourth, fifth, and sixth grid points of the discretized space. I know via broadcasting I can obtain the absolute valued squared data by abs.sol

I tried surface(range(-5,5), sol.t, abs.(sol)), but this resulted in an error.

A surface is either defined explicitly as z=f(x,y), (x,y)\in D\subset\mathbb{R}^2 or parametric: (u,v)\to(x(u,v), y(u,v), z(u,v)), with (u,v) in a a rectangular region or a disk, etc.
Your sol has length(sol.t)= length(sol.u)=19, and to each sol.t[j], j=1:19,
corresponds a vector sol.u[j] of length 6.

If you associate to each sol.t[j] the absolut values,
abs.(sol.u[j])=(abs(sol.u[j][1]), abs(sol.u[j][2]), ... abs(sol.u[j][6])
you get a map s\in \mathbb{R}\to\mathbb{R}^6, and it doesn’t define a surface, but a curve in a 6-dimensional space.

The spatial data, (abs(sol.u[j][1]), abs(sol.u[j][2]), ... abs(sol.u[j][6]) are slices of the surface at time sol.u[j].

Well, if you take:

using Plots
plotlyjs()
 x=sol.t
 y = range(-5, 5, length=6)
 z= [abs(sol.u[i][j])  for j=1:6, i=1:19]
 surface(x, y, z)

you’ll get the surface.
Since you didn’t give any information on the six y-values for each sol.t[i], I defined them as being equally distanced.

Just in case you weren’t aware of the existing integrations with Plots.jl: Plot Functions