**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.