\subsection{Some test runs and results}
For the modifications
\begin{equation}
  f(x,y)={\frac {\sqrt {x}\left (1-x\right ){y}^{1/10}\left (1-y\right
  )^ {2}}{3}}  
\end{equation}
and  $N=10$ made to the listed code, we have the following relative
errors:
For $f$,
$$
\left [\begin {array}{ccccc}  0.03479& 0.02597& 0.02024& 0.02052&
 0.01858\\\noalign{\medskip}- 0.007716&- 0.003380&- 0.007902&- 0.003937
&- 0.006782\\\noalign{\medskip}- 0.01327&- 0.008304&- 0.01731&- 0.02282
&- 0.02138\\\noalign{\medskip} 0.01171& 0.004684& 0.005910& 0.0001523&-
 0.0004227\\\noalign{\medskip} 0.001112&- 0.009258&- 0.007214&- 0.01015
&- 0.01108\end {array}\right ]
$$
\noindent
for $\partial f/\partial x$:
$$
\left [\begin {array}{ccccc}  0.02671&- 0.008716& 0.02433& 0.009778&-
 0.01909\\\noalign{\medskip}- 0.03521&- 0.07032& 0.03061&- 0.006469&-
 0.07817\\\noalign{\medskip} 0.004837&- 0.04405&- 0.06126&- 0.04526&
 0.05540\\\noalign{\medskip} 0.02411& 0.001753&- 0.01669&- 0.04427&
 0.04858\\\noalign{\medskip}- 0.003292&- 0.01564&- 0.01042&- 0.04122&
 0.01453\end {array}\right ]
$$
\noindent
for $\partial f/\partial y$:
$$
\left [\begin {array}{ccccc}  0.07000& 0.2547& 0.07463&- 0.2004&-
 0.1303\\\noalign{\medskip} 4.226& 2.744& 4.528& 4.940& 4.182
\\\noalign{\medskip}- 0.6513&- 0.4244&- 0.6238&- 0.5545&- 0.5229
\\\noalign{\medskip}- 0.03197& 0.05759&- 0.01336&- 0.04826&- 0.03390
\\\noalign{\medskip} 0.1819& 0.1774& 0.2088& 0.1813& 0.1747\end {array}
\right ]
$$
\noindent
on the grid with x-values
$$
\left [\begin {array}{ccccc}  0.01& 0.05750& 0.1050& 0.1525& 0.2000
\\\noalign{\medskip} 0.01& 0.05750& 0.1050& 0.1525& 0.2000
\\\noalign{\medskip} 0.01& 0.05750& 0.1050& 0.1525& 0.2000
\\\noalign{\medskip} 0.01& 0.05750& 0.1050& 0.1525& 0.2000
\\\noalign{\medskip} 0.01& 0.05750& 0.1050& 0.1525& 0.2000\end {array}
\right ]
$$
\noindent
and y-values
$$
\left [\begin {array}{ccccc}  0.01& 0.01& 0.01& 0.01& 0.01
\\\noalign{\medskip} 0.05750& 0.05750& 0.05750& 0.05750& 0.05750
\\\noalign{\medskip} 0.1050& 0.1050& 0.1050& 0.1050& 0.1050
\\\noalign{\medskip} 0.1525& 0.1525& 0.1525& 0.1525& 0.1525
\\\noalign{\medskip} 0.2000& 0.2000& 0.2000& 0.2000& 0.2000\end {array}
\right ]
$$

With the ``nicer'' function
\begin{equation}
f(x,y)={\frac {\sqrt {x}\left (1-x\right ){y}^{3/5}\left (1-y
\right )^{2}}{3}}  
\end{equation}
and $N=10$, using the same $(x,y)$ grid, we have the relative errors
$$
\left [\begin {array}{ccccc}  0.0147& 0.00644&0& 0.00147&0
\\\noalign{\medskip}- 0.00757&0&- 0.00646&0&- 0.00524
\\\noalign{\medskip}- 0.0102&- 0.00641&- 0.0150&- 0.0175&-
 0.0202\\\noalign{\medskip}- 0.00261&- 0.0114&- 0.00445&-
 0.0156&- 0.0144\\\noalign{\medskip} 0.00249&- 0.00545&-
 0.00425&- 0.00744&- 0.0103\end {array}\right ]
$$

$$
\left [\begin {array}{ccccc} - 0.0100&- 0.0281& 0.00918&0&-
 0.0380\\\noalign{\medskip}- 0.0309&- 0.0695& 0.0319&- 0.00540
&- 0.0754\\\noalign{\medskip} 0.00896&- 0.0420&- 0.0589&-
 0.0438& 0.0583\\\noalign{\medskip} 0.00532&- 0.00748&- 0.0293
&- 0.0520& 0.0375\\\noalign{\medskip}- 0.00254&- 0.0143&-
 0.00932&- 0.0408& 0.0138\end {array}\right ]
$$
\noindent and
$$
\left [\begin {array}{ccccc} 0& 0.0259&- 0.00289&- 0.0379&-
 0.0304\\\noalign{\medskip}- 0.0569& 0.00499&- 0.0700&- 0.0954
&- 0.0693\\\noalign{\medskip} 0.0378&- 0.0295& 0.0273&0&-
 0.00349\\\noalign{\medskip} 0.0994&- 0.0181& 0.0622& 0.0867&
 0.0665\\\noalign{\medskip} 0.102& 0.0185&- 0.0764&- 0.00744&-
 0.00688\end {array}\right ]
$$
%
for $f$, $\partial f/\partial x$ and $\partial f/\partial y$, respectively.


With the non-separable function
 \begin{equation}
f(x,y)=\sqrt [4]{{x}^{2}+{y}^{2}}   
 \end{equation}
and $N=10$, again on the same grid, we get the following results.
Notice that this function is {\sl not} in $L_{\alpha}$ on the unit
square, so poor results can be expected.
 \begin{equation}
\left [\begin {array}{ccccc} - 0.0841&0& 0.00616&- 0.0153&
 0.00894\\\noalign{\medskip}- 0.0124&0&- 0.00289&- 0.0198&-
 0.00438\\\noalign{\medskip}- 0.00616&- 0.00289&- 0.0104&-
 0.0372&- 0.0168\\\noalign{\medskip}- 0.0179&- 0.0248&- 0.0209
&- 0.0345&- 0.0179\\\noalign{\medskip} 0.0134&- 0.00658&-
 0.00421&- 0.0160&- 0.00752\end {array}\right ]
 \end{equation}
%
 \begin{equation}
\left [\begin {array}{ccccc} - 0.0875& 0.324&- 0.235& 0.0705&
 0.254\\\noalign{\medskip}- 1.80& 0.231&- 0.205& 0.0173& 0.161
\\\noalign{\medskip}- 0.103& 0.388&- 0.396&- 0.0199& 0.271
\\\noalign{\medskip} 2.29& 0.538&- 0.317&- 0.00657& 0.274
\\\noalign{\medskip}- 4.27& 0.367&- 0.362&- 0.0546& 0.265
\end {array}\right ]
 \end{equation}
%
 \begin{equation}
\left [\begin {array}{ccccc} - 0.138&- 0.110&- 2.13&- 13.6&-
 16.1\\\noalign{\medskip} 0.289& 0.377& 0.147&- 0.357&- 0.499
\\\noalign{\medskip}- 0.215&- 0.268&- 0.210&- 0.162&- 0.200
\\\noalign{\medskip} 0.118& 0.0173& 0.0606& 0.130& 0.101
\\\noalign{\medskip} 0.302& 0.207& 0.170& 0.179& 0.158
\end {array}\right ]
 \end{equation}
%
for $f$, $\partial f/\partial x$ and $\partial f/\partial y$, respectively.


Next, an $L_{\alpha}$ function with a radial singularity at $(0,0)$.
 \begin{equation}
f(x,y)=-{y}^{3/2}\left (-1+y\right )\left (x-1\right )+\sqrt 
[4]{{x}^{2}+{y}^{2}}\left (1-x\right )\left (1-y\right )y
 \end{equation}
With $N=10$ and the above grid, we have the relative errors of 
 \begin{equation}
\left [\begin {array}{ccccc}  0.0162&0&- 0.00502&- 0.00410&-
 0.00363\\\noalign{\medskip} 0.0613&0&- 0.00777&- 0.00398&-
 0.00640\\\noalign{\medskip}- 0.189&- 0.0103&- 0.0136&-
 0.0236&- 0.0176\\\noalign{\medskip}- 0.339&- 0.00622&-
 0.00869&- 0.0210&- 0.00871\\\noalign{\medskip}- 0.271&-
 0.00739&- 0.00498&- 0.0149&- 0.00923\end {array}\right ]
 \end{equation}
%
 \begin{equation}
\left [\begin {array}{ccccc}  0.00691&- 0.0170&- 0.00881&
 0.0128&- 0.00371\\\noalign{\medskip} 0.00528&- 0.0246&
 0.00359&- 0.00451&- 0.0265\\\noalign{\medskip}- 0.0814&
 0.0151&- 0.0434&- 0.0257& 0.0143\\\noalign{\medskip}- 0.141&
 0.0427&- 0.0436&- 0.0203& 0.0177\\\noalign{\medskip}- 0.122&
 0.0339&- 0.0305&- 0.0136& 0.0126\end {array}\right ]
 \label{eqn:test3-dx}
 \end{equation}
%
 \begin{equation}
\left [\begin {array}{ccccc}  0.00662& 0.0194&- 0.0130&-
 0.0349&- 0.0257\\\noalign{\medskip} 0.190&- 0.0975&- 0.108&-
 0.0931&- 0.0527\\\noalign{\medskip} 0.785& 0.0236&- 0.0421&-
 1.31&- 0.0367\\\noalign{\medskip}- 0.214&- 0.00650&- 0.0553&
- 0.120&- 0.211\\\noalign{\medskip}- 1.09&- 0.0171&- 0.0172&-
 0.0546&- 0.0402\end {array}\right ]
 \label{eqn:test3-dy}
 \end{equation}
%
for $f$, $\partial f/\partial x$ and $\partial f/\partial y$, respectively.
The (4,1) and (5,1) entries in (\ref{eqn:test3-dx}) correspond to
small exact values;  similarly, the poor entries of
(\ref{eqn:test3-dy}) also correspond to small exact answers.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
As yet another nice test, consider
 \begin{equation}
f(x,y)=\sqrt {x}\sqrt [3]{1-x}{y}^{1/10}\left (1-y\right 
)^{2}\sin(2\,\pi \,x){e^{y}}
 \end{equation}
with $N=10$ and the rest as coded.  For the function interpolation, we
have the relative errors
 \begin{equation}
\left [\begin {array}{ccccc}  0.0385& 0.0251& 0.0168&
 0.0212& 0.0203\\\noalign{\medskip}- 0.00452&- 0.00680&-
 0.00738&- 0.00935&- 0.00358\\\noalign{\medskip}- 0.0360&
- 0.00851&- 0.0147&- 0.0233&- 0.0178\\\noalign{\medskip}-
 0.0115& 0.00174&0&- 0.00476&0\\\noalign{\medskip}-
 0.0240&- 0.0127&- 0.00784&- 0.00993&- 0.0114\end {array}
\right ]
 \end{equation}
Along the diagonal from (0,0) to (1,1), the relative error is similar,
with a maximum of 5\%.

Without the random error, we obtain
 \begin{equation}
\left [\begin {array}{ccccc}  0.0667& 0.0289& 0.0251&
 0.0317& 0.0324\\\noalign{\medskip} 0.0362&- 0.00341&-
 0.00738&0&0\\\noalign{\medskip} 0.0360&- 0.00340&0&0&0
\\\noalign{\medskip} 0.0508& 0.0104& 0.00753& 0.0143&
 0.0146\\\noalign{\medskip} 0.0337&- 0.00543&- 0.00784&0&0
\end {array}\right ]
 \end{equation}
The results along the diagonal are again comparable, with a maximum
relative error of 5\%.
By increasing $N$ to 13, we obtain the relative error of 
 \begin{equation}
\left [\begin {array}{ccccc}  0.00770& 0.0193& 0.00836&
 0.0212& 0.0203\\\noalign{\medskip}- 0.0203&- 0.0119&-
 0.0148&- 0.0141&- 0.00716\\\noalign{\medskip}0& 0.0119&
 0.00735& 0.00465& 0.0107\\\noalign{\medskip}- 0.0161&-
 0.00522&- 0.0151&- 0.00952&- 0.00364\\\noalign{\medskip}
- 0.0168&- 0.00907&- 0.00784&- 0.00497&- 0.00380
\end {array}\right ]
 \end{equation}
Along the diagonal, the maximum relative error is still 5\%.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


If we define $\rho $ by 
 \begin{equation}
\rho = r^{\frac{1}{2}}\sin (\theta) = \frac{y}{(x^{2}+y^{2})^{1/4}}
 \end{equation}
with $r$ the radial distance from the origin, and $\theta $ measured
from the x-axis as usual, then the function
 \begin{equation}
f(x,y) = \rho (1-x)(1-y) - \sqrt{y} (1-x)(1-y)
 \end{equation}
is in $L_{\alpha}$ on the unit square.
With $N=10$ and the above grid, we have the relative errors of 
 \begin{equation}
\left [\begin {array}{ccccc}  0.0449& 0.0110& 0.00163&- 0.00160
&0\\\noalign{\medskip}- 0.00602&- 0.00295&- 0.00644&- 0.00514&-
 0.00583\\\noalign{\medskip}- 0.196&- 0.0115&- 0.0145&- 0.0231&
- 0.0203\\\noalign{\medskip} 0.723& 0.0216& 0.00365&- 0.0134&-
 0.0119\\\noalign{\medskip}- 0.357&- 0.0106&- 0.00529&- 0.0152&
- 0.0110\end {array}\right ]
 \end{equation}
%
 \begin{equation}
\left [\begin {array}{ccccc} 0&0&- 0.138&- 0.0471&- 0.0336
\\\noalign{\medskip}- 0.0185&- 0.0148&- 0.00225&0&- 0.0356
\\\noalign{\medskip}- 0.0546&0&- 0.0380&- 0.0265& 0.0276
\\\noalign{\medskip}- 0.0129& 0.0317&- 0.0394&- 0.0269& 0.0198
\\\noalign{\medskip}- 0.118& 0.0282&- 0.0277&- 0.0139& 0.0164
\end {array}\right ]
 \end{equation}
%
 \begin{equation}
\left [\begin {array}{ccccc} - 0.0231& 0.107&- 0.00605&- 0.0703
&- 0.0516\\\noalign{\medskip} 0.784& 0.0404& 0.108& 0.228&
 0.443\\\noalign{\medskip}- 1.64&- 0.0688&- 0.0650&- 0.0425&-
 0.0436\\\noalign{\medskip} 0.576& 0.0197&- 0.0133&- 0.0341&-
 0.0394\\\noalign{\medskip} 4.28& 0.140& 0.0575& 0.00482&0
\end {array}\right ]
 \end{equation}
%
for $f$, $\partial f/\partial x$ and $\partial f/\partial y$, respectively.
The 0.443 entry in the y-partial corresponds to a small entry.


For the same $f$ and $N$, but without introducing the random error,
the relative errors of 
 \begin{equation}
\left [\begin {array}{ccccc}  0.0449& 0.0146& 0.0114& 0.0112&
 0.0114\\\noalign{\medskip}- 0.0181&0&0&0&0\\\noalign{\medskip}
- 0.0911&- 0.00576&- 0.00242&0&- 0.00135\\\noalign{\medskip}
 0.974& 0.0402& 0.0146& 0.00672& 0.00683\\\noalign{\medskip}-
 0.163&- 0.00453&0&- 0.00305&- 0.00220\end {array}\right ]
 \end{equation}
%
 \begin{equation}
\left [\begin {array}{ccccc}  0.00489&0&- 0.0199& 0.0637&
 0.0122\\\noalign{\medskip}- 0.00308& 0.00269&- 0.00225&0&0
\\\noalign{\medskip}- 0.00781&- 0.00404& 0.00422&0&0
\\\noalign{\medskip} 0.0703& 0.00633&- 0.00262&0&0
\\\noalign{\medskip}- 0.0273& 0.00471&- 0.00346&0& 0.00410
 \end {array}\right ]
 \end{equation}
%
 \begin{equation}
\left [\begin {array}{ccccc} - 0.00771& 0.0273& 0.0121& 0.00502
& 0.00939\\\noalign{\medskip} 0.809& 0.0827& 0.0841& 0.136&
 0.352\\\noalign{\medskip}- 1.81&- 0.101&- 0.0542&- 0.0504&-
 0.0591\\\noalign{\medskip} 0.734& 0.0344& 0.0133& 0.0137&
 0.00657\\\noalign{\medskip} 4.79& 0.199& 0.0719& 0.0434&
 0.0410\end {array}\right ]
 \end{equation}
%
are obtained for $f$, $\partial f/\partial x$ and $\partial f/\partial
y$, respectively.  Notice the large errors for small $x$, in the first
columns of the matrices.
Along the diagonal, for $x=y=0.1:0.1:0.9$, the largest relative error
is 0.5\%.  
% 	left boundary error:
For $x=0.02, y=0.1:0.1:0.9$ (near the left boundary),  we
get 
$$
[ 0.065, 0.024,
0.511,-1.32, 0.00302,
3.32,-3.32,- 0.815,- 8.07]
$$  
for the values of the relative error.  
%	top boundary error
For $x=0.1:0.1:0.9, y=0.98$ (the top boundary), the r.e. is 
$$
[ 1.29, 11.7,- 1.13,- 0.473,- 0.267,- 0.179,-0.137,- 0.107,- 0.0893] 
$$
while for $y=0.1:0.1:0.9, x=0.98$ (the right boundary) it is much improved:
$$
[0,0, 0.00283,0,0, 0.00579,0,- 0.00137,- 0.00598]
$$
%
%	bottom boundary error
For $x=0.1:0.1:0.9, y=0.02$
(the bottom boundary), we also have the much better error
$$
[ 0.00896, 0.00779,
 0.00723, 0.00665,
 0.00658, 0.00659,
 0.00633, 0.00634,
 0.00622]
$$

For the above $f, N,$ without random error and with $h =
\frac{\pi}{\sqrt{N}}$ we have the following errors near the boundaries:
On the top,
 \begin{equation}
[- 0.0108,- 0.151,- 0.0946,- 0.0339, 0.000233,
 0.0123, 0.00746,- 0.00496, 0.00899]
 \end{equation}
on the bottom,
 \begin{equation}
[ 0.00197,- 0.000588, 0.000203, 0.000743,
 0.000122,- 0.000729,- 0.000223, 0.00203,-
 0.00152]
 \end{equation}
on the left, 
 \begin{equation}
[-0.0453, 0.170,-
 0.0583,- 0.39,-
 0.0761, 0.828,
 0.556,- 2.692,
 3.44]
 \end{equation}
and on the right
 \begin{equation}
[ 0.0102,- 0.00506,- 0.0203,- 0.00692,- 0.00441
,- 0.0237,- 0.0172,- 0.0321,- 0.0210]
 \end{equation}
The graphs comparing the exact answer and the interpolant are shown in
the figures.