Orthogonal polynomial interpolation for 2 dimension.

   [f, esys, eeval, edom] = interporth2(x, Fxt, domain, errorsurf)

   f = interporth2(x, Fxt) interpolates an expression by using the
   Chebyshev or Legendre points; the approximation is
   P(x, t) = sum(i = 0:n-1)sum(j = 0:n-1)a_{i, j}Px_i(x)Pt_j(t), where
         A = [a_{0, 0} a_{0, 1} ... a{0, N}]
             [a_{1, 0} a_{1, 1} ... a{1, N}]
             [           ...            ]
             [a_{N, 0} a_{N, 1} ... a{N, N}], where N = n-1. 

Inputs (required)
   x         = independent tau variable(itau object).
   Fxt       = two variable function F(x, t) (char).     

Inputs (optional)
   domain    = rectangular domain (double vector: [x1 x2 y1 y2]).
   errorsurf = bolean variable to surf the error.

   f         = interpolation coeficients (double matrix).
   esys      = error in the linear system (double vector).
   eeval     = maximum error in the evaluation (double scalar).
   edom      = error in the full domain (double matrix).

See also
  interporth, interporthinc and interporth2inc.