determined the real numbers $c_{10},\cdots c_{(k-1)0}$ we obtain the real number $c_{k0}$ by solving $$\label{eikatk}
d_{k0}^2(x_k,y)=\frac{1}{4}\sum_{ij} a_{ij}(x_k) d^ 2_{k0,x_i}(x_k,y)d^ 2_{k0,x_j}(x_k,y).$$ for $c_{k0}$. Continuing this procedure for $N$ interpolation points we get a polynomial of the form $$\label{dN1}
d_{N0}^2(x,y)=d^ 2_{A^ {-1}(y)}(x,y)+\sum_{j=1}^Nc_{j0}\Pi_{r=0}^j\Pi_{l=1}^n(x_l-x^r_l)^2d^ 2_{A^ {-1}(x^j)}(x,y).$$ with $N$ real numbers $c_{j0}$ obtained recursively by plugging in $d^2_{j0}$ with one degree of freedom $c_{j0}$ into .

Analogous constructions are done to approximate $x\rightarrow d^2(x,x^j)$ for $k=1,\cdots ,N$ with $$\label{dNk}
d_{Nk}^2(x,x^k)=d^ 2_{A^ {-1}(y)}(x,x^k)+\sum_{j=1}^Nc_{jk}\Pi_{r=0}^j\Pi_{l=1}^n(x_l-x^r_l)^2d^ 2_{A^ {-1}(x^j)}(x,x^k),$$ with $c_{jk}$ computed analogously. The construction of the functions $d^2_{N0},\cdots , d^2_{NN}$ suffices to approximate $d^2$ (we do not need to synthesize these functions into one function, for example by a Lagrangian polynomial interpolation). Note that for $j=0,\cdots N$ the function $d^2_{Nk}$ satisfies the equation $$\label{bdpk}
\begin{array}{ll}
d^2(x,x^k)=\frac{1}{4}\sum_{ij}a_{ij}(x)d^2_{x_i}(x,x^k)d^2_{x_j}(x,x^k)\\
\\
\mbox{ with boundary condition }\\
\\
d^2(x,x^k)=0~\mbox{ iff }~x=x^k.  
\end{array}$$ at all interpolation points $x^0,\cdots x^N$ by construction.

Note that in the preceding construction no restrictions on the choice of the interpolation points are made. This does not mean that one may search for an optimal choice of interpolation points and improve efficiency and convergence. We are free to choose a certain set of interpolation points (for example Chebyshev nodes). But these 