Looking for a way to powerfully descritize a PDE (Julia)

The problem I'm trying to solve is this PDE,

I wrote the buck of the solver as follows (mostly based on the Euler Explicit):

for n in 2:N-1
    wl[n]= w[n] + Δx*v[n]
    ϕl[n]= ϕ[n] + Δx*ρ[n] 
    ρl[n] = ρ[n] - (Δt/Δx)*(dif_nt(ρ,n)*v[n]) + (Δt)*w[n]*ρ[n]
    vl[n] = v[n] - (Δt/Δx)*(v[n]*w[n]) + (μ*Δt/(ρ[n]*Δx))*(dif_nt(w,n)) + (c₀^2*Δt/ρ[n]*Δx)*(dif_nt(ρ,n)) + (Δt/τ)*(V(ρ[n])-v[n])
end

In which the dif_nt is simply the forward distance u[n+1] - u[n] for a variable. nt stands for no-time.

function dif_nt(v, n)
    return v[n+1] - v[n]
end

How could I use a library to find a good approximation to (pl[n],vl[n],ϕ[n],wl[n]), e.i., these are the next-step values in time?