Plots.jl not showing negative values in 3D

I am solving a 3D differential equation. My output data is in the form of a NxNxN dense Matrix for each timestep. With values ranging from -1 to 1. I plan on animating the plot trough time so that leaves 3 dimensions.

The standard surface/contour/plot/heatmap from Plots.jl all produce the same image. They dont allow me to change the color limits and dont display negative numbers.
Heres a MWE of the plotting.

function initial(N)
    as = LinRange(0,pi,N)
    u0 = [sin(x)^2*sin(y)^2*sin(z)^2 for x in as, y in as, z in as]
    return u0
end
u0 = initial(50)
plot(u0)#this works okay 

plotu0

plot(-u0)#this doesnt work

plot_4
The plot technically works for positive values but it is difficult to understand what you see. And the plot doesnt display negative values at all.
What would be the best way to fix this?
And in general is there a better way to plot 3D data?

Such volume rendering is only available for Plots.jl gr() backend and from the GR docs I believe it is using a Maximum Intensity Projection (MIP) technique - tbc.

So when negating the data example provided, the 0s will be the largest amplitudes projected.

Plotting 2D slices of the 3D volume may be easier to read.