I am trying to plot the solution vector `sol`

from the solver of an `ODEProblem`

which is a 1432 element vector of 31x1 vectors. To do this, I attempted

```
plotlyjs()
y= sol.t
z= [abs(sol.u[i][j])^2 for j=1:N₁, i=1:size(sol)[2]]
surface(x, y, z)
```

where `N₁=31`

However, this never finishes compiling. What is the proper way of plotting my solution? Below is the full MWE.

```
using DifferentialEquations, LinearAlgebra, Plots, SparseArrays
N₁=31 # Number of waveguides
γ=-1 # Nonlinear term strength parameter
h=0.6 # 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]*Tridiagonal(fill(1,N₁-1), fill(-2,N₁), fill(1,N₁-1))*u + γ*u.^3)
u0 = sech.(x);
u0 = map(Complex,u0)
tspan = (0.0,100)
prob = ODEProblem(g₁,u0,tspan, [h])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
plotlyjs()
y= sol.t
z= [abs(sol.u[i][j])^2 for j=1:N₁, i=1:size(sol)[2]]
surface(x, y, z)
```