The immediate issue is that ModuleApply, which is what is called by TestModule(3), returns not a module but the result of a computation.  Modify the procedure to return thismodule.  Also, there is a bug in CodeTools:-Profiling:-PrintProfiles; you can work around it by replacing that line by

printf("%s",debugopts('procdump' = run));

Using modules as dynamic objects is possible, however, the newer Maple construct, an object (implemented as a specialization of a module) is generally better.

Change

ModelSet(1):= Model1;

to

ModelSet(1):= eval(Model1);

See the help page last_name_eval for a discussion of why this is necessary.

A simple way to see understand what happens during parallel transport on a constant latitude of a sphere is to consider parallel transport on the cone that is tangent to the sphere at that latitude. A cone is flat, so the change in angle is equal to the angle of the cut from a flat sheet that is used to form the cone.

One approach is to return the output as a string, then use StringTools:-RegSubs to remove the lhs.

printf("%s", StringTools:-RegSubs("^[^=]*= "="", C(x^2+y^2+z^2,'output=string'))):
x * x + y * y + z * z;

Write V as a vector (not a list). LinearAlgebra:-VectorNorm can be used to compute the norm of the vector.

V := <vx,vy>:
-1/2*C*LinearAlgebra:-VectorNorm(V,2)*V;

Here is one approach

myproc := proc(A,B,C)
local L;
   L := [A,B,C];
   if L :: 'list(algebraic)' then
      # option 1
   elif L :: listlist then
      if numelems(A) = 2 then
          # option 2
      elif numelems(A) = 3 then
          # option 3
      end if;
   elif L :: 'list(Vector)' then
      if numelems(A) = 2 then
          # option 4 (may want to check others for consistency)
      else
          # option 5 (ditto)
      end if;
    else
        error "unknown types of arguments"
    end if;
end proc;

Consider using mtaylor rather than two calls to series.

mtaylor(f(x,y),[x,y],3);

One could use the optional third argument in a call to collect to apply a function to each coefficient of a polynomial.

Directly assigning values, or in this case, lists of values, to variables that are used in the expressions of interest is not generally the best strategy with Maple. More effective is to assign a set of equations that define the numeric values of the variables and then use them to evaluate the expressons. For example

params := [ k = 1, m = 2 ]:
y := k/m^2:
eval(y, params);

Your situation is more complicated in that you have multiple values for each parameter. I assume that you want to evaluate the derived expressions with all parameters taking, say, the second value in their respective lists. You could accomplish that by doing something like

params := [ k = [1,2], m = [2,3] ]:
getparams := proc(k) local eq; [seq(lhs(eq) = op(k, rhs(eq)), eq = params)]; end proc:
eval(y, getparams(2));

You could do

foldl((a,b)->mods(a*b,m),x,y,z);
           mods(mods(x y, m) z, m)

The simple answer is yes, the profiling feature doesn't work with gotos.

That aside, you really should not be using gotos, and not just for Dijkstra's structural concerns.  Maple gotos are for the most part undocumented.  The kernel implementation is not efficient (a goto statement is implemented via a search through the statements for the label, put more unexecuted statements between them and it will take longer).  The argument to the goto is evaluated and the target cannot be a local variable, so there are namespace issues.

Followup

The namespace issues can be avoided by using strings as labels. A little experimentation indicates that the search for the label appears to progress by first checking the preceding statement, then the following statement, then the statement before the preceding statement, etc.

The content of that library does not include squarefree; it appears to only have Auxiliary.  How did you create the respository?

Use forward quotes.  For good measure you should also use them around units:

with(Units[Standard]):
with(ScientificConstants):
c := evalf(Constant('c', 'units'));

Thanks for reporting this.  I'm submitting a bug report.

Use the factors command.  For example

   y := (a+b)^3:
   f := factors(y);
                   f := [1, [[a + b, 3]]]
   map(L -> `$`(op(L)), f[2]);
                    [a + b, a + b, a + b]

See help on @@:  (f@@3)(x);