  We study quartic matrix models with partition function Z[E,J]=\int dM
\exp(trace(JM-EM^2-(\lambda/4)M^4)). The integral is over the space of
Hermitean NxN-matrices, the external matrix E encodes the dynamics, \lambda>0
is a scalar coupling constant and the matrix J is used to generate correlation
functions. For E not a multiple of the identity matrix, we prove a universal
algebraic recursion formula which gives all higher correlation functions in
terms of the 2-point function and the distinct eigenvalues of E. The 2-point
function itself satisfies a closed non-linear equation which must be solved
case by case for given E. These results imply that if the 2-point function of a
quartic matrix model is renormalisable by mass and wavefunction
renormalisation, then the entire model is renormalisable and has vanishing
\beta-function.
  As main application we prove that Euclidean \phi^4-quantum field theory on
four-dimensional Moyal space with harmonic propagation, taken at its
self-duality point and in the infinite volume limit, is exactly solvable and
non-trivial. This model is a quartic matrix model, where E has for N->\infty
the same spectrum as the Laplace operator in 4 dimensions. Using the theory of
singular integral equations of Carleman type we compute (for N->\infty and
after renormalisation of E,\lambda) the free energy density
(1/volume)\log(Z[E,J]/Z[E,0]) exactly in terms of the solution of a non-linear
integral equation. Existence of a solution is proved via the Schauder fixed
point theorem.
  The derivation of the non-linear integral equation relies on an assumption
which we verified numerically for coupling constants 0<\lambda\leq (1/\pi).
