For an N'th order truncated multivariable Taylor series expansion for f(x,y) about (x₀,y₀), 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:

temp.mw

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 ?

I want to explore multivariable function approximations using truncated Taylor series.

Mathematically, for a function f(x,y) and using operator notation for the partial derivatives, where e.g. D_x² f(x,y) denotes the second partial derivative of f wrt x evaluated at (x,y), we can write the N'th order truncated Taylor series for f around (x₀,y₀) as

I want to make a Maple-function for this expression, and try

P := (x,y,x0,y0,N) -> sum(1/factorial(n)*sum(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):

where f(x,y) is a previosly defined Maple-function.

My P function fails, and the reason why it fails is related to the "D" operator in the "sum".

Please take a look at the following code-snippet:

Output (9) is as expected, but output (8) is not !!
I would expect output (8) to be equal to the sum of output (9), i.e. to be equal to (-1/2).

Please illuminate why I don't get the sum of the sequence (9) as my output (8).

Is there any alternative to Mathematica's "RegionFunction" in Maple, with applications of controlling the plotting domain for surface plots ? This question has been asked previosly, and given a workaround by "Markiyan Hirnyk", suggesting to leave the function that is to be plotted undefined outside of the plot region. Something like:

f1 := proc(x,y) if x^2 + y^2 > 0 and x^2 + y^2 < 4 then 8 - x^2 - y^2 end if end proc:
f2 := proc(x,y) if x^2 + y^2 > 0 and x^2 + y^2 < 4 then x^2 + y^2 end if end proc:
plot3d({f1,f2}, -2..2, -2..2) 

The result is clearly not good looking. You can increase the number of points, but we need a better way to do this.

BTW, the similar question was asked in 2010. Are there any better alternatives, ~ 10 years later ?

Simple example to illustrate the desired functionality:
Say we have a 2D vector function which describes the position of a particle

r := t -> <5*cos(Pi*t), 5*sin(Pi*t)>

We want to define the velocity and acceleration as functions, so we could do something like

v := t -> <diff(r(t)[1], t), diff(r(t)[2], t)>

The problem now is that we cannot call our velocity function with numeric arguments.
A simple solution is to call the function via "subs", as in

subs(t = 2, v(t))

but IMHO, this is not very elegant and I guess inefficient. Is there a command that enables for pulling out the evaluated result from diff such that it can be used directly as a functional expression ? I.e., I want to be able to call

v(2)

directly, without having to do substitutions.

EDIT:

I found that you can do

v := <diff(r(t)[1], t, diff(r(t)[2], t)>
v := unapply(v, t)

but please provide your recommendations. Thanks

I have started to use Maple to test my calculations for a complex variable course.

The example is a complex integrand, and integration about an ellipse centered at origin.
Integrand has function 

f(z) = z*exp(a*z) / (z²+1)²

C:  z(t) = cos(t) + i*2*sin(t)

where a is a real constant. My idea is shown in the attatched workbook, i.e. 

dz = -sin(t) + i*2*cos(t) dt
g(t) = (-sin(t) + i*2*sin(t)) * f(cos(t) + i*sin(t))
int(g(t), t = 0..2*Pi)

PS, I've done similar code previously, but only for circle contours at center (see attachement)

The analytical answer is pi*a*sin(a), and Maple gives this correct for circular contours, but not for the ellipse.

So, what I'm asking is not how to solve the problem, but how to perform (arbitrary) contour integrals in Maple.

attachForMaplePrime.mw