Heat trace for Laplace type operators with non-scalar symbols

For an elliptic selfadjoint operator P=−[uμν∂μ∂ν+vν∂ν+w] acting on a fiber bundle over a compact Riemannian manifold, where uμν,vμ,w are N×N-matrices, we develop a method to compute the heat-trace coefficients ar which allows to get them by a pure computational machinery. It is exemplified in any even dimension by the value of a1 written both in terms of uμν=gμνu,vμ,w or diffeomorphic and gauge invariants. We also address the question: when is it possible to get explicit formulae for ar?