The fluid here has been assumed as single component perfect gas i.e. it obeys the equation $p=ρ R T$, the thermal conductivity is assumed as a constant, so the equation set is: NDSolve[{D[ρg[t, x], t] + D[ρg[t, x] u[t, x], x] == 0, ρg[t, x] D[u[t, x], t] + ρg[t, x] u[t, x] D[u[t, x],x] == -D[ρg[t, x] te[t, x], x], ρg[t, x] D[te[t, x], t] + ρg[t, x] u[t, x] D[te[t, x],x] == -ρg[t, x] te[t, x] D[u[t, x], x] + D[te[t, x], x, x], te[0, x] == 298, te[t, -0.5] == 298, te[t, 0.5] == 298, ρg[0, x] == (1 - x^2), u[0, x] == 0, u[t, -0.5] == 0, u[t, 0.5] == 0}, {ρg[t, x], te[t, x], u[t, x]}, {t, 0, 1}, {x, -0.5, 0.5}] After I ran the code, I only got the warning message `NDSolve::ndsz` and `NDSolve::eerr`, I've checked the the equations for times and I think they are correct, and the initial and boundary conditions are also simple and seems to be reasonable (at least from the perspective of physics). So…What's wrong with it?…Well, to tell you the truth, what I really want to ask is, does `NDSolve` lack the ability to solve system of partial differential equations? Oh, someone may feel strange that there's no boundary condition for `ρg[t, x]`, that's because, I found that only four boundary conditions are necessary for the solving of the equations though I don't know the exact reason (I found it in times of trial when I set `ρg[0, x]` as a constant ).