then quite similar as in the stage for multiindices of order $3$. Therefore we give a very short description. Starting with the multiindex $\delta^0$ there is a multiindex $\beta^k$ of order $m$ (i.e. $|\beta^k|=m$) such that $\delta^0-\beta^k=1_i$ for some index $i$. Then we get successive approximations $$d^2_{\delta^0 k}(x,y)=d^2_{M(\beta_m)}(x,y)+\sum_{r=0}^{k}c_{\delta^0}^r\Pi_{l=0}^{r-1}(x-x^l)^{\delta^0+{\bf 1}}\frac{1}{\delta^0!}(x-x^r)^{\delta^0},$$ where the real numbers $c_{\delta^0}^k$ are succesively determined by plugging in the function $x\rightarrow d^2_{\delta^0 k}(x,y)$ into the equation $$\partial^ {\beta^k}d^2=\frac{1}{4}\partial^ {\beta^k}\left( \sum_{lm}a_{lm}(x)d^2_{x_l}d^2_{x_m}\right) ,$$ evaluated at the interpolation point $x^k$ (Note that $\partial^{\beta^k}=\partial^ {\delta^0-1_i}$). After $N+1$ steps we get the approximation function $x\rightarrow d^2_{\delta^0 N}(x,y)$. Having defined $x\rightarrow d^2_{\delta^l N}(x,y)$ for $l=0,\cdots p-1$ the next multiindex $\delta^r$ may be such that there is an multiindex $\beta^h$ of order $m$ such that $\delta^{r}-\beta^h=1_i$ for some index $i$. We may then define $x\rightarrow d^2_{\delta^p k}(x,y)$ $$d^2_{\delta^p k}(x,y)=d^2_{M(\beta_m)}(x,y)+\sum_{r=0}^{k}c_{\delta^p}^r\Pi_{l=0}^{r-1}(x-x^l)^{\delta^p+{\bf 1}}\frac{1}{\delta^p!}(x-x^r)^{\delta^p},$$ and determine the constants $c_{\delta^p}^r$ by plugging in the function $x\rightarrow d^2_{\delta^r k}(x,y)$ into the equation $$\partial^ {\beta^h}d^2=\frac{1}{4}\partial^ {\beta^h}\left( \sum_{lm}a_{lm}(x)d^2_{x_l}d^2_{x_m}\right) ,$$ and evaluate at $x^k$. Finally, we get the approximation of order $m+1$, namely $$\label{betam}
d^2_{M(\beta_{m+1})}=d^2_{\delta^{R_{m+1}} N}(x,y).$$ Note that this approximation satisfies the eikonal equation and all its derivatives up to order $m+1$, i.e. all equations $$\partial^{\alpha}_x d^2=\frac{1}{4}\partial^{\alpha}_x \left( \sum_{lm}a_{lm}(x)d^2_{x_l}d^2_{x_m}\right)$$ with $|\alpha|\leq m+1$ at all 