I have a stiff system error, any solution I have tried tries to evaluate it, but then stops(don't know why). P[x_, y_, t_] = e[x, y, t]/(γ - 1) ; e[x_, y_, t_] = (γ - 1) ρ[x, y, t]/(μ mu ); kb T[x, y, t]; η = 10^-6; cp = 5/2 kb/(μ mu); Rgas = 8.3144598; cv = 5/2 kb/(μ mu) - Rgas; γ = cp/cv; g = 28.02*9.81; μ = 0.6163328197226503`; mu = 1.66053904*10^-27; kb = 1.38064852*10^-23; sol = NDSolve[{D[T[x, y, t], y] == 0.41 g ρ [x, y, t] T[x, y, t]/P[x, y, t], D[P[x, y, t], y] == g ρ[x, y, t], T[x, 0, t] == 5770, ρ[x, 0, t] == 1.42*10^-7*1.408*10^3}, {ρ, T}, {y, 0, 4*10^6}, {x, 0, 12*10^6}, {t, 0, 10}]; sol1 = NDSolve[{ D[ρ[x, y, t]*u[x, y, t], t] == -D[ρ[x, y, t]*u[x, y, t]*u[x, y, t] + P[x, y, t], x] - D[ρ[x, y, t]*u[x, y, t]* v[x, y, t] - η (D[u[x, y, t], y] + D[v[x, y, t], x]), y], D[ρ[x, y, t]*v[x, y, t], t] == -D[ρ[x, y, t]*v[x, y, t]* u[x, y, t] - η (D[v[x, y, t], x] + D[u[x, y, t], x]), y] - D[ρ[x, y, t]*v[x, y, t]*v[x, y, t] + P[x, y, t], y] + g ρ[x, y, t], D[ρ[x, y, t], t] == -D[ρ[x, y, t]*u[x, y, t], x] - D[ρ[x, y, t]*v[x, y, t], y], D[e[x, y, t], t] == -D[u[x, y, t]*e[x, y, t], x] - D[v[x, y, t]*e[x, y, t], y] - P[x, y, t]*(D[u[x, y, t], x] - D[v[x, y, t], y]), v[0, y, t] == v[12*10^6, y, t], u[0, y, t] == u[12*10^6, y, t], T[0, y, t] == T[12*10^6, y, t], ρ[0, y, t] == ρ[12*10^6, y, t], e[x, 0, t] == 3.83767261162, v[x, 4000000, t] == 0, v[x, 0, t] == 0, (D[u[x, y, t], y] /. y -> 0) == 0, (D[u[x, y, t], y] /. y -> 4000000) == 0, v[x, y, 0] == 0, u[x, y, 0] == 0, T[x, y, 0] == 5770 + 0.00835414960707927 y, ρ[x, y, 0] == 1.42*10^-7*1.408*10^3 + 7.3561137493644*10^-10 y }, {u, v, T, ρ}, {x, 0, 12000000}, {y, 0, 4000000}, {t, 0, 100}] The error is as follows: NDSolve::ndsz: At t == 0.0005560756736763107`, step size is effectively zero; singularity or stiff system suspected. I also get this error: NDSolve::eerr: Warning: scaled local spatial error estimate of 332.45447642223894` at t = 0.0005560756736763107` in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 13 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options. I do not know if the first error causes the second or visa-versa, my area is rather large, but I do not know how to fix the gridspace problem, as I said, when I try to solve it with possible solutions(increasing maximum stepsize for example) the thing stops for no given reason, perhaps this has something to do with ram or maximum running time. As you can see, I'm not that experienced with solving stuff in higher than 1 dimension, I hope this is a resolvable problem. EDIT: When I add a MaxStepFraction of 1/100 the code runs, but the cell keeps running and doesn't produce results. But unlike my comments earlier it keeps running and doesn't stop randomly