Two dimensional Taylor series

                 jn   kn  hx^j  hy^k  d(j)  d(k) 
 F(x+hx,y+hy) = sum  sum  ----  ----  ----- ----- F(x,y)
                j=0  k=0   j!    k!   dx(j) dy(k) 

         d(j)  d(k) 
 denote  ----- ----- F(x,y) as  Fxxxyy   for j=3, k=2 
         dx(j) dy(k) 

 
The first few terms are: (remember hy^0 = 1, 0! = 1)
 F(x+hx,y+hy) = F                        ! j=0, k=0  thru nj=0, nk=0

              + hx Fx                    ! j=1, k=0
              + hy Fy                    ! j=0, k=1
              + hx hy Fxy                ! j=1, k=1  thru nj=1, nk=1

              + hx^2/2! Fxx              ! j=2, k=0
              + hx^2/2! hy/1!  Fxxy      ! j=2, k=1
              + hx^2/2! hy^2/2! Fxxyy    ! j=2, k=2
              + hy^2/2! Fyy              ! j=0, k=2
              + hx/1! hy^2/2! Fxyy       ! j=1, k=2  thru nj=2, nk=2

              + hx^3/3! Fxxx             ! j=3, k=0
              + hx^3/3! hy/1! Fxxxy      ! j=3, k=1
              + hx^3/3! hy^2/2! Fxxxyy   ! j=3, k=2
              + hx^3/3! hy^3/3! Fxxxyyy  ! j=3, k=3
              + hy^3/3! Fyyy             ! j=0, k=3
              + hx/1! hy^3/3! Fxyyy      ! j=1, k=3
              + hx^2/2! hy^3/3! Fxxyyy   ! j=2, k=3  thru nj=3, nk=3



 The above can be placed in a two dimensional matrix, A :

          j=0        j=1         j=2        j=3        j=nj
 k=0    1/(0! 0!)  1/(1! 0!)   1/(2! 0!)  1/(3! 0!)
 k=1    1/(0! 1!)  1/(1! 1!)   1/(2! 1!)  1/(3! 1!)
 k=2    1/(0! 2!)  1/(1! 2!)   1/(2! 2!)  1/(3! 2!)
 k=3    1/(0! 3!)  1/(1! 3!)   1/(2! 3!)  1/(3! 3!)

 k=nk                                                 1/(nj! nk!)

 two dimensional matrix, B :

         j=0        j=1         j=2        j=3       
 k=0      F          Fx          Fxx        Fxxx
 k=1      Fy         Fxy         Fxxy       Fxxxy
 k=2      Fyy        Fxyy        Fxxyy      Fxxxyy
 k=3      Fyyy       Fxyyy       Fxxyyy     Fxxxyyy 

 two dimensional matrix, C :

         j=0        j=1         j=2        j=3       
 k=0   hx^0 hy^0   hx^1 hy^0   hx^2 hy^0   hx^3 hy^0
 k=1   hx^0 hy^1   hx^1 hy^1   hx^2 hy^1   hx^3 hy^1
 k=2   hx^0 hy^2   hx^1 hy^2   hx^2 hy^2   hx^3 hy^2
 k=3   hx^0 hy^3   hx^1 hy^3   hx^2 hy^3   hx^3 hy^3

 Then the sum of the element by element product, A*B*C, is the Taylor series
 expansion of F(x+hx,y+hy) in terms of F(x,y) and its derivatives

 F(x,y) is a given function that can be evaluated numerically.

 Fx(x,y), Fxx(x,y), Fxxx(x,y), ... can be computed numerically as one variable
                                   see  nderiv.c  and  nderiv.out
                                   and use F(x+a*hx,y) type terms

 Fy(x,y), Fyy(x,y), Fyyy(x,y), ... can be computed numerically as one variable
                                   see  nderiv.c  and  nderiv.out
                                   and use F(x,y+b*hy) type terms

 Fxy(x,y) and all mixed derivatives require numerical partial differentiation

 for example Fxy(x,y) = 1/(4hx hy) ( F(x+hx,y+hy) - F(x+hy,y-hy) -
                                     F(x-hx,y+hy) + F(x-hx,y-hy) )
                                     + error O(hx hy)

 for higher order, for example Fxxyy(x,y) =
   1/(hx^2 hy^2) (F(x+hx,y+hy) + F(x-hx,y+hy) + F(x+hx,y-hy) + F(x-hx,y-hy) +
                  F(x,y) - 2(F(x+hx,y) + F(x-hx,y) + F(x,y+hy) + F(x,y-hy)) )