I need to solve a system of DE's in 2 dimensions over time, The problem is that mathematica complains when I have a boundary condition of the same order as the differential order of the DE's, so I tried using the finiteelement method since it worked for other people who had similar problems, but it just returns itself as if nothing ever happened P[x_, y_, t_] = e[x, y, t]/(γ - 1) ; e[x_, y_, t_] = (γ - 1) ρ[x, y, t]/(μ mu ) kb T[x, y, t]; 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; 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], y], D[ρ[x, y, t]*v[x, y, t], t] == -D[ρ[x, y, t]*v[x, y, t]* u[x, y, t], 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], (D[u[x, y, t], x] /. x -> 0) == (D[u[x, y, t], x] /. x -> 12000000), (D[v[x, y, t], x] /. x -> 0) == (D[v[x, y, t], x] /. x -> 12000000), D[u[0, y, t], y] == D[u[12000000, y, t], y], D[v[0, y, t], y] == D[v[12000000, y, t], y], D[u[0, y, t], t] == D[u[12000000, y, t], t], D[v[0, y, t], t] == D[v[12000000, y, t], t], (D[T[x, y, t], x] /. x -> 0) == (D[T[x, y, t], x] /. x -> 12000000), (D[ρ[x, y, t], x] /. x -> 0) == (D[ρ[x, y, t], x] /. x -> 12000000), D[T[0, y, t], y] == D[T[12000000, y, t], y], D[ρ[0, y, t], y] == D[ρ[12000000, y, t], y], D[T[0, y, t], t] == D[T[12000000, y, t], t], D[ρ[0, y, t], t] == D[ρ[12000000, y, t], t], (D[e[x, y, t], x] /. x -> 0) == (D[e[x, y, t], x] /. x -> 12000000), D[e[0, y, t], y] == D[e[12000000, y, t], y], D[e[0, y, t], t] == D[e[12000000, y, t], 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}, Method -> "FiniteElement" ]