What about the following solution?

        remove(has,rhs(aa),omega^(-sigma+1))
+factor(select(has,rhs(aa),omega^(-sigma+1)));

The following short example, using one of your functions, may (ask otherwise) get you started on how to write the module itself:

M := module()
   option package;
   export Quadrance;
   Quadrance := (a1,a2) -> (a1[1] - a2[1])^2 + (a1[2] - a2[2])^2;
end module:
M:-Quadrance([1,2],[2,4]);   # An example of use

To be able to load that module with the with-command, save it with savelib(M), where I assume that the variable savelibname equals the library path you give; see the help page module[package].

After having loaded the module with the with-command, the long form M:-, as used in the example above, is no longer needed; see the help page UsingPackages.

Perhaps the following function is useful:

f := (expr) -> mtaylor(expr,select(x -> type(x,indexed) and op(0,x) = B,indets(expr)),2):

Applied to any one of your polynomials, it will linearize it with respect to any indexed variable with the name B.

Perhaps the following is useful:

removeProducts := proc(x)
   local expr := expand(x);
   if type(expr,Or(`&*`(
      anything,
      Or(identical(a),identical(a)^posint),
      Or(identical(b),identical(b)^posint)
   ),`&*`(
      Or(identical(a),identical(a)^posint),
      Or(identical(b),identical(b)^posint)
   ))) then 0
   elif type(expr,Not(atomic)) then map(removeProducts,expr)
   else expr
   end if
end proc:
removeProducts(f0);
removeProducts(f1);
removeProducts(f2);

PS: The type specification can quite possibly be improved upon.

Perhaps I am misunderstanding what you want to do, but what about the following?

createFunc := (f::list,vars::list) -> unapply(f,vars[]):

Example of use:

F := createFunc([y,y*z-x,-15*x*y-x*z-x],[x,y,z]):
F(a,b,c);
F(1,1,1);

I assume that you by "not solving for the 'D' operator" mean "not performing the substitution for the D's". Perhaps the following is helpful:

eval(convert(junk,diff),x[1](T[0],T[1],T[2]) = R(T[1],T[2])*sin(omega*T[0]+phi(T[1],T[2])));

The following recipe works for me:

Expand all sections.
Go to the very bottom of the worksheet.
Insert a new section.
Mark from the start of the new section to the very bottom of the worksheet, and press delete.

Is the following useful?

expr := diff(lambda(t-t1),t);
expr := subs(lambda = sin,expr);
convert(expr,diff);

Or perhaps simpler, using eval rather than subs:

expr := diff(lambda(t-t1),t);
eval(expr,lambda = sin);

The line vec := map(f,vec) in the approx procedure is illegal, as vec is a parameter of the procedure; it cannot be assigned anything from within the procedure.

If your problem is due to overly long program lines, then try to break them (at proper places, of course) using Shift+Enter (line break).

A plot of a complex function, i.e., a function from C to C, requires four real dimensions, so unless you have special mental capabilities :-), it is impossible to visualize.

Perhaps the help page on DEtools[symgen] would be helpful.

Use L := seq(i,i=1..4), not L = seq(i,i=1..4).

Warning (note added): The answer is misconceived, see my comment 'Probably mistaken' below.

I am surprised to find that there seems to be no command in the LieAlgebras package concerning isomorphisms. As a lesser substitute, I have come up with the following solution (at least I think it is):

with(DifferentialGeometry):
with(LieAlgebras):
L1 := LieAlgebraData([[e2, e3] = e3, [e2, e4] = -e4, [e3, e4] = -e1],[e1,e2,e3,e4],alg1);
L2 := LieAlgebraData([[e2, e3] = e1, [e4, e2] =  e2, [e3, e4] = -e3],[e1,e2,e3,e4],alg2);
DGsetup(L1);
DGsetup(L2);
A := Matrix(4,symbol = a);
Query(alg1,alg2,A,indets(A),"Homomorphism")[4];

According to Example 2 on the help page for Query, the above construction finds all homomorphisms that sends alg1 to alg2. As all resulting matrices are singular, there would seem to be no isomorphisms between the two algebras.

Perhaps I am missing the point, but the sum of odd integers from 1 to 2*k + 1 equals s(k) = (k + 1)^2, i.e., a square. There exists no integer k for which s(k) = 3375; the sums closest to 3375 are s(57) = 3363 and s(58) = 3480, respectively.