Suppose I use the DifferentialEquations.jl package and find a solution structure sol.

Is there a simple way to find the derivative of the various variables in sol? In other words: suppose I have 6 states in an ODE model, and want to plot the time derivative of state nr. 4?

I guess your differential equations are of the form dx_i/dt = f_i(x,t). Therefore, you can obtain the time derivative directly by evaluating the differential equations.

You just add a second argument with the order of derivative you want:

julia> begin
using DifferentialEquations
f(u,p,t) = 1.1*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-10, abstol=1e-10)
end;
julia> sol(1.0) # solution at t = 1.0
1.50208301197696
julia> sol(1.0, Val{1}) ) # First derivative of the solution at t = 1.0
1.6522913131746269
julia> 1.1 * sol(1.0) ≈ sol(1.0, Val{1}) # We know for this differential equation, that this should match the derivative for any t
true