For an N'th order truncated multivariable Taylor series expansion for f(x,y) about (x0,y0), you can write
f := (x,y) -> 1/(2 + x*y^2):
P := (x,y,x0,y0,N) -> add(1/factorial(n)*add(binomial(n,k)*D[1$(n-k),2$k](f)(x0,y0)*
(x-x0)^(n-k)*(y-y0)^k, k=0..n), n=0..N):
The procedure works well, but it is slow for large orders.
Comparing with Maple's mtaylor function:
The built-in mtaylor is ~1400 times faster at N=30, and increasing as N gets larger.
What is the easiest way to speed up my P procedure ?