taylor2.c running 
poly2_print nj=4, nk=4 
26.000000  24.000000  21.000000  17.000000  
25.000000  22.000000  18.000000  13.000000  
23.000000  19.000000  14.000000   8.000000  
20.000000  15.000000   9.000000   2.000000  

trim 
poly2_print nj=4, nk=4 
26.000000  24.000000  21.000000  17.000000  
25.000000  22.000000  18.000000   0.000000  
23.000000  19.000000   0.000000   0.000000  
20.000000   0.000000   0.000000   0.000000  

poly2_eval P1 at x=1, y=1, is 215 
poly2_evaln P1 at x=1, y=1, is 215 
poly2_eval P1 at x+hx=1.25, y+hy=1.5, is 409.281 
poly2_evaln P1 at x+hx=1.25, y+hy=1.5, is 409.281 
derivative with respect to x = 209, using pderiv 
derivative with respect to x = 209, is P2 
poly2_print nj=3, nk=4 
25.000000  22.000000  18.000000   0.000000  
46.000000  38.000000   0.000000   0.000000  
60.000000   0.000000   0.000000   0.000000  

derivative with respect to y = 194, using pderiv 
derivative with respect to y = 194, is P3
poly2_print nj=4, nk=3 
24.000000  42.000000  51.000000  
22.000000  36.000000   0.000000  
19.000000   0.000000   0.000000  
 0.000000   0.000000   0.000000  

2nd derivative with respect to x = 204, using pderiv 
2nd derivative with respect to x = 204, is P4
poly2_print nj=2, nk=4 
46.000000  38.000000   0.000000   0.000000  
120.000000   0.000000   0.000000   0.000000  

2nd derivative with respect to y = 180, using pderiv 
2nd derivative with respect to y = 180, is P5
poly2_print nj=4, nk=2 
42.000000  102.000000  
36.000000   0.000000  
 0.000000   0.000000  
 0.000000   0.000000  

derivative with respect to x and y = 96, using pderiv 
derivative with respect to x and y = 96, is P6
poly2_print nj=3, nk=3 
22.000000  36.000000   0.000000  
38.000000   0.000000   0.000000  
 0.000000   0.000000   0.000000  

3rd derivative with respect to x = 120, using pderiv 
3rd derivative with respect to x = 120, is P7
poly2_print nj=1, nk=4 
120.000000   0.000000   0.000000   0.000000  

3rd derivative with respect to y = 102, using pderiv 
3rd derivative with respect to y = 102, is P8
poly2_print nj=4, nk=1 
102.000000  
 0.000000  
 0.000000  
 0.000000  

derivative with respect to xx and y = 38, using pderiv 
derivative with respect to xx and y = 38, is P9
poly2_print nj=2, nk=3 
38.000000   0.000000   0.000000  
 0.000000   0.000000   0.000000  

derivative with respect to x and yy = 36, using pderiv 
derivative with respect to x and yy = 36, is P10
poly2_print nj=3, nk=2 
36.000000   0.000000  
 0.000000   0.000000  
 0.000000   0.000000  

all fourth derivative matrices should be zero 
build test poly2 about (x,y) 
just factorials A 
poly2_print nj=4, nk=4 
 1.000000   1.000000   0.500000   0.166667  
 1.000000   1.000000   0.500000   0.166667  
 0.500000   0.500000   0.250000   0.083333  
 0.166667   0.166667   0.083333   0.027778  

factorials and hx, hy  AB 
poly2_print nj=4, nk=4 
 1.000000   0.500000   0.125000   0.020833  
 0.250000   0.125000   0.031250   0.005208  
 0.031250   0.015625   0.003906   0.000651  
 0.002604   0.001302   0.000326   0.000054  

build F manually, automatic in taylor2_eval 
F(x,y), Fx(x,y) ... in  F 
poly2_print nj=4, nk=4 
215.000000  194.000000  180.000000  102.000000  
209.000000  96.000000  36.000000   0.000000  
204.000000  38.000000   0.000000   0.000000  
120.000000   0.000000   0.000000   0.000000  

full matrix ABF 
poly2_print nj=4, nk=4 
215.000000  97.000000  22.500000   2.125000  
52.250000  12.000000   1.125000   0.000000  
 6.375000   0.593750   0.000000   0.000000  
 0.312500   0.000000   0.000000   0.000000  

ABF(x,y,xh,yh) = 409.281 same as poly2_eval(x+xh,y+yh) = 409.281 

taylor2_eval(x,y,xh,yh) = 409.281 same as poly2_eval(x+xh,y+yh) = 409.281

taylor2.c finished