because $$\begin{array}{ll}
d^ 2(x,y)=&\frac{1}{4}\sum_{ij}a_{ij}d^ 2_{x_i}d^ 2_{x_j}=\frac{1}{4}\sum_{ij}\frac{1}{2}\left( a_{ij}+a_{ji}\right) d^ 2_{x_i}d^ 2_{x_j}\\
\\
&+\frac{1}{4}\sum_{ij}\frac{1}{2}\left( a_{ij}-a_{ji}\right) d^ 2_{x_i}d^ 2_{x_j}=\frac{1}{4}\sum_{ij}\frac{1}{2}\left( a_{ij}+a_{ji}\right) d^ 2_{x_i}d^ 2_{x_j},
\end{array}$$ so we can always substitute the matrix $a_{ij}$ by its symmetrization $\frac{1}{2}\left( a_{ij}+a_{ji}\right)$ without affecting the solution $d^ 2$.

In [@Ka] we have seen that for $C^{\infty}$ coefficient functions $x\rightarrow a_{ij}(x)$ and $x\rightarrow b_i(x)$ and if some boundedness conditions of the derivatives are satisfied the fundamental solution has the pointwise valid form $$\label{WKBrep}
p(t,x,y)=\frac{1}{\sqrt{2\pi t}^n}\exp\left(-\frac{d^2(x,y)}{2 t}+\sum_{k= 0}^{\infty}c_k(x,y)t^k\right),$$ where the functions $x\rightarrow c_k(x,y),~k\geq 0$ are solutions of recursively defined linear first order equations for each $y$. These equations can be solved by methods of characteristics or approximated by regular polynomial interpolation methods outlined in [@Ka2]. In the computation of the WKB-coefficients $d^2$ and $c_k,~k\geq 0$ the recursive relations for $c_{k+1}$ involve second order derivatives of $c_k$, and therefore implicitly derivatives of order $2k$ of the squared metric $d^2$. Hence it is of great interest to compute not only $d^2$ but also its derivatives up to a given order with high accuracy. The present work shows how his can be accomplished. In Section 2 we recapture some facts about the connection of the geodesic equation , and equation , and prove 