LU decomposition is numerically stable, for the most part. You can run into problems with floating point arithmetic if you're not careful. This just has to do with the way you solve the problem. By breaking up the matrix into an upper and lower triangular part and solving the resulting triangular system, you often avoid things like powers or square roots and just reduce everything to basic addition and multiplication which tend to have more stable numerical properties.
As for speed, well in this small case of N=4 there is probably very little speed increase since LU is O(N^{3}), although this can be improved depending on symmetries of the matrix. Maybe for N=5 I might write out the exact formulas but beyond that I would just a LU solver since the amount of code for an LU solver is fairly compact.
But, as you correctly point out, if you are never changing the matrix A then you are probably safe with writing it out this way.
I just come from a numerical fluids background so seeing exact formulas makes me uneasy. And in my work since the number of grid points N tends to be variable so general formulas are not available. Plus I've noticed a lot of programmers tend to not know numerical algorithms. In fairness, I don't blame them. All the numerical computations courses I've taken and TAed are very boring and don't make you do anything fun and numerical analysis has this reputation of being dry. If you actually made them write stuff like a Navier-Stokes solver, you would get them interested.
By the way if you want to learn more about numerical linear algebra, which is the cornerstone of most scientific and high performance computing, I personally enjoy Trefethen and Bau's book. Although it's aimed at a mathematical audience and assumes such.