When plotting **DifferentialEquations** solutions, the *plotdensity* keyword does not seem to work together with *line_z*. See code below where data interpolation within **Plots** does not take place, when we color the solution by the simulation time array (`line_z=sol0.t`

):

```
using DifferentialEquations, ParameterizedFunctions, Plots
lorenz = @ode_def begin
dx = σ*(y - x); dy = x*(ρ - z) - y; dz = x*y - β*z
end σ ρ β
u0=[1.,1.,1.]; tspan=(0.,60.); p=[10.,28.,8/3];
prob = ODEProblem(lorenz, u0, tspan, p)
sol0 = solve(prob, reltol=1e-9)
# plot(sol0, vars = (1, 2, 3)) # display interpolation is on but line_z not working with this syntax
plot((sol0[1,:],sol0[2,:],sol0[3,:]), line_z=sol0.t, plotdensity=10_000, color=:redsblues)
```

Is this possible to do using Plots?