@Mac Dude 

mseries works in your example, but as mentioned in the answer,  x=0  cannot be replaced with x.
In your example you must use:

map(vvmseries,TM,[x=0,xp=0,y=0,yp=0,dl=0,dp=0,dkQ1=0],3);

On the other side, iterated series does not always work;
An example:

restart;
(* Test case for mtaylor *)

vvmseries:=proc(f,X::list(name=anything),n::posint)
local t;
  eval( f, `=`~(lhs~(X), rhs~(X)+t*~(lhs~(X)-rhs~(X))) );
  eval(convert(series(%, t, n), polynom),t=1);
end proc:
mymseries:=proc(f,vars,norder)
      local i,fs:=f;
          for i from 1 to numelems(vars) do
              fs:=convert(series(fs,vars[-i],norder),polynom);
          end do;
          return fs;
end proc:

f:=(x+y)^x:
f1:=vvmseries(f,[x=0,y=0],5):
f2:=mymseries(f,[x,y],5):
eval([f,f1,f2],[x=0.0001, y=0.001]);

[.9993189875,  .9993189874,  .7786567398]

Best regards,

V.A.

[edited example]

@one man 

1. The only problem is that the thread is about parametrization of a surface and you insist in posting curves (animated or not).
Your curves are not bad, but their place is not here.

2. I do not have a general method. In your example:
(x1^2+x2^2-0.4)^2+(x3+sin(x1*x2+x3))^4-0.1=0;

a) start with a parametrization U(t), V(t) of the curve U^2 + V^4 - 0.1 = 0.
b) define g(v,w) as the unique root z of    z + sin(w + z) = v
     (In Maple it us a RootOf)
c) The parametrization is
     x := sqrt(U(t) + 0.4) * cos(u)
     y := sqrt(U(t) + 0.4) * sin(u)
     z := g(V(t), x*y).

Here also implicitplot3d is your main ingredient to see the surface.
Compare with a real surface parametrization (exact or approx).

@Markiyan Hirnyk 

1. It is similar and simpler. So, have you some series in mind for this?

3. Am I supposed to write a procedure for generalized series in several variables from scratch?

@Markiyan Hirnyk 

For series (in a single variable) of a non-analytic function Maple uses "generalized series". Such series do not exist for several variables (or I don't know about them) so this is all that can be obtained.
E.g. what should be the series for (x-y)^(1/3) at x=0,y=0?

@Markiyan Hirnyk 

It works for me. Try a restart.

@Carl Love 

Ok, but it's still a bug. The function does not have discontinuities, so the option should be superfluous.
The antiderivative has a removable singularity at 0 so for Maple it remains to execute the subtraction (FTOC).
Note that f:=diff(sin(x^2)/x,x); also has a removable singularity at x=0 but int works ok.

@Markiyan Hirnyk 

Compare:

mseries(x^x*sin(y), [x=0,y=1], 3);
and
mtaylor(x^x*sin(y), [x=0,y=1], 3);

@Christopher2222 
Of course it works but the bug refers to the symbolic computation.

A curve (having 1 parameter) cannot substitute a surface parametrization which has 2 parameters.
I will not continue, I understood.

@one man 

It seems that you found a very serious BUG.

h:=foo(x) / erf(x)+a;
int(h,x);
      a*x

This also happens if erf(x) is replaced by erf(x^(1/3))  etc.
For a=0 the integral remains unevaluated.

@one man 

If a parametrization is there, then please use it to plot the surface, instead of using implicitplot3d. As I said, it would be very useful but unfortunately I don't see it and I don't think it will be easy to obtain.

I think you ask too much from Maple, the PDF for sin(X) seems to be hard to compute.

Let's take a much simpler distribution where Maple also fails but which can be easily computed.

X := RandomVariable(Uniform(0, 4*Pi));

Maple says that the PDF for sin(X) is

f := piecewise(t < 0, 0, t = 0, 5/(4*Pi), t < 1, 1/(sqrt(-t^2+1)*Pi), 1 <= t, 0)

but the correct PDF is 2*f.

@one man 

You have parametrized two curves (red an green, lying on the surface), but you have not parametrized the surface.

@one man 

No need for other examples of animated curves, you already posted dozens.
If you can, show us a real (approximate) parametrization for this surface, this would be useful.