We show that the average characteristic polynomial P_n(z) = E [\det(zI-M)] of
the random Hermitian matrix ensemble Z_n^{-1} \exp(-Tr(V(M)-AM))dM is
characterized by multiple orthogonality conditions that depend on the
eigenvalues of the external source A. For each eigenvalue a_j of A, there is a
weight and P_n has n_j orthogonality conditions with respect to this weight, if
n_j is the multiplicity of a_j. The eigenvalue correlation functions have
determinantal form, as shown by Zinn-Justin. Here we give a different
expression for the kernel. We derive a Christoffel-Darboux formula in case A
has two distinct eigenvalues, which leads to a compact formula in terms of a
Riemann-Hilbert problem that is satisfied by multiple orthogonal polynomials.