multiindex $\beta^1$ from the list of multiindices of order $3$ (i.e. $|\beta^1|=3$) where we may assume that $\beta^1-\gamma^1=1_k$ for some multiindex $\gamma^1$ with $|\gamma^1|=2$ and some index $k$. An analogous construction as in the case of $\beta^0$ can be done. The only difference is that we start with $d^2_{\beta^0 N}(x,y)$ instead of $d^2_{N0}(x,y)$. We get an approximation of the form $$d^2_{\beta^1 N}(x,y)=d^2_{\beta^0 k}(x,y)+\sum_{k=0}^Nc_{\beta^1}^k\Pi_{l=0}^{k-1}(x-x^l)^{\beta^1+1}\frac{1}{\beta^1!}(x-x^k)^{\beta^1}.$$ where the real numbers are computed recursively by plugging the current approximation into the equation $$\label{beta0}
\begin{array}{ll}
\partial^{(\beta^1-\gamma^1)}_xd^2=\partial^{(\beta^1-\gamma^1)}_x\left(\frac{1}{4}\sum_{ij}a_{ij}(x)\frac{\partial d^2}{\partial x_i}\frac{\partial d^2}{\partial x_j} \right),
\end{array}$$ evaluating at the current interpolation point and solving for the currently undetermined real number $c_{\beta^1}^k$. Doing this for all the multiindices of order $3$ in the list above we get the approximation $$d^2_{M(\beta_3)}(x,y):=d^2_{\beta^N N}(x,y).$$ Note that by construction the function $x\rightarrow d^2_{M(\beta_3)}(x,y)$ satisfies the equation and all its first order derivative equations $$\partial^{i}_xd^2=\frac{1}{4}\partial^{i}_x\left( \sum_{lm}a_{lm}(x)d^2_{x_l}d^2_{x_m}\right) ,~~1\leq i\leq n,$$ at all interpolation points $x^0=y,x^1,\cdots, x^N$. This completes the stage of construction for multiindices of order $3$. Next assume that the construction for the approximation $$x\rightarrow d^2_{M(\beta_m)}(x,y)$$ of order $m$ has been completed. Then we may list the multiindices of order $m+1$, i.e. consider a list of multiindices $\delta^0,\delta^1,\cdots, \delta^{R_{m+1}}$ such that $|\delta|=m+1$. The procedure is 