- You must give numerical values to parameters, e.g.

omega:=1; M:=1; N:=1; N1:=1; Delta:=1;

(you probably did).

- Choose first a smaller interval, e.g.

plots[odeplot](res,[[t,(Re(w(t)))]],0..5,axes=boxed,titlefont=[SYMBOL,14],font=[1,1,18],color=black,linestyle=1,tickmarks=[3, 4],font=[1,1,14],thickness=2,titlefont=[SYMBOL,12]);

(for larger intervals it will take several minutes)

It works.

I think that unapply does not work properly if its first argument is a procedure.

So,

g := unapply(x -> theta*x, theta);

also does not work.
unapply is not builtin, so it can be examined.

A curious fact is that f and g are apparently the same, according to print or lprint,
but they have distinct internal structures according to
dismantle(eval(f)); dismantle(eval(g));

Why not simply use Diff instead of diff? Then op() works.

Otherwise, here is another (non-recursive) function:

undiff:=proc(f)
local u:=f, X:=NULL;
while op(0,u)=diff do X:=op(2,u),X; u:=op(1,u) od;
u,[X]
end;

The expressions being so complicated, Maple has problems recognizing some zero coefficients.

Workaround:

u := evecs[1][1];
simplify(numer(u)):
u1:=simplify(series(%,varepsilon=0)):
simplify(denom(u)):
u2:=simplify(series(%,varepsilon=0)):
evalc(series(u1/u2,varepsilon=0));

You have symbolic expressions containing floats. You should simplify them first do avoid possible catastrophic cancelations.

Use e.g.

Ez := simplify(Ez);

and compare with the unsimplified version.

Let G be the list of the polynomials obtained by Basis. E.g.

G:=[y^3+x*y, x*y^2+x^2, x^2*y-2*y^2+x, x^3-3*x*y]:

mon:=proc(p) local t; coeffs(p,indets(p),t); t end:
mon ~ (G);

convert(%,set); # remove duplicates

createModule2 := proc(A::Matrix(square))
    local dim;
    dim := RowDimension(A);
    module()
        export det;
        det := (x::Matrix(1..dim,1..dim)) -> Determinant(x);
        det :=subs('dim'=dim, eval(det));   ###
    end module
end proc;

Try with

Digits:=30;

Your polynomial ferrai contains 3 symmetric sums as coefficients. By knowing these, you can't hope to recover all 4 fundamental symmetric sums (which are in fact c1,...,c4)  simply because 3<4.

BTW, if you have a symmetric sum s, you can express it in terms of c1,...,c4 using

simplify(s, [x1+x2+x3+x4=c1,...,x1*x2*x3*x4=c4]);

pochammer function appears only because the result must be also valid for i=0.

(For integer i>0, pochhammer(1-i, i)=0.)

How would you answer the following scalar version?

Given two numbers a and b, calculate the number c in terms of a and b.

Example for your 2nd example:

s:=sum((sum(f[t](x[i,k],a[k]),k=1..n)-y[i, t])*(sum(f[t](x[i, k], a[k]),k=1..n)-y[i,t]),i=1..m);

# diff(s, a[p]) = ?

s11:=op([1,1],s); # identify the sum containing a[p]
s11new:=subs([(k=1..n) = (k=1..p-1)], s11) +
                subs([(k=1..n) = (k=p+1..n)], s11) +
                subs([(k=1..n) = (k=p..p)], s11):
snew:=subs(s11=s11new,s):
diff(snew,a[p]);

I have attached a worksheet.
An element of the algebra is represented as Vector(v1,...,vn,v0); also included a display procedure (optional)

e[i] &x e[j]  = 0 if not defined otherwise.

I have modified f because the indices in your definition exceed n; you also have provided two incompatible variants.

algebra-eiej-code-1.mw  (edit: updated)

You want Groebner[Basis], not LinearAlgebra[Basis] ! So, use:

G := Groebner[Basis](K, 'tord', degrevlex(r,u,v));

It is not very clear what are you trying to obtain.

Mathematically a polynomial map is a vector function (a,b,c) |--> (r,u,v) having polynomial components.

In Maple you may represent it by a list [r,u,v] or vector Vector([r,u,v]) provided that the indeterminates (i.e. variables) are a,b,c.

Or, you may consider the procedures:
p:=unapply( [r,u,v], a,b,c );
or
p:=unapply( Vector([r,u,v]), a,b,c );