@Axel Vogt Thank you!
You forgot "return" before 'procname'.
I still wonder what goes on in plot and plot3d.

@Axel Vogt There is a difference in the way plot and plot3d handles modp:
restart;
modp(5.4321,2); #An error as expected
plot(x->modp(x,2),-3..3); #Works, surprisingly
plot('modp(x,2)',x=-3..3); #Works, surprisingly
plot3d((x,y)->modp(x,2),-3..3,-1..1); #Doesn't work
plot3d('modp(x,2)',x=-3..3,y=-1..1); #Works

@Axel Vogt There is a difference in the way plot and plot3d handles modp:
restart;
modp(5.4321,2); #An error as expected
plot(x->modp(x,2),-3..3); #Works, surprisingly
plot('modp(x,2)',x=-3..3); #Works, surprisingly
plot3d((x,y)->modp(x,2),-3..3,-1..1); #Doesn't work
plot3d('modp(x,2)',x=-3..3,y=-1..1); #Works

Your solution works, but I wonder why since
F(5.78);
results in the error

Error, (in F) invalid argument for modp or mods

as does also
modp(5.78,2);

Your solution works, but I wonder why since
F(5.78);
results in the error

Error, (in F) invalid argument for modp or mods

as does also
modp(5.78,2);

I fail to understand why you expect answers of type numeric in what you are doing. In fact I don't understand what you are trying to do.
Your equations eq1 and eq2 could be obtained in a somewhat more readable form like this:

restart;
egf := exp(2*x*z-z^2);
test1r := subs(x=dummy,z=x,egf);
test1r2 := solve(test1r=f,x);
b := Matrix([[a1,a2],[a3,a4]]);
invb := b^(-1);
testing1:=<f|f^2>.invb.Matrix([[x,1],[1,x^2]]).b; #Row vector
eq1 := test1r2[1] = testing1[1];
eq2 := test1r2[1] = testing1[2];

To understand why you want to linearize it would be good to see a full example. Is expression supposed to mean expression = 0 and is it part of a system involving a derivative of g?

Since diff(f(r), r$2) is already linear in f you only need to linearize the rest.
The command
mtaylor(exp(-f)*(1 - g)^2,[f,g],2);
linearizes about zero, so it is the same as:
mtaylor(exp(-f)*(1 - g)^2,[f=0,g=0],2);

When making the decision back in 1993 (I think it was) the price and licensing certainly favored Maple. So did the open code. The first version we got was Maple V, Release 2, but we started using it in our teaching with Release 3. That students could have Maple on their own computer at home was essential.

@CarlitosVillalbaGalea You can do
solO1:=subs(sol3,1-C1(t)-C2(t)-C3(t));
allvalues(solO1); #Length of output exceeds limit of 1000000
#But it seems to be OK by checking with
eval(%,params); #Assuming that params has been defined already

@CarlitosVillalbaGalea You can do
solO1:=subs(sol3,1-C1(t)-C2(t)-C3(t));
allvalues(solO1); #Length of output exceeds limit of 1000000
#But it seems to be OK by checking with
eval(%,params); #Assuming that params has been defined already

I replaced all floating point numbers by 3 or -3 depending on their signs. As an example, if the number 1.23456*10^12 appeared in your system it would have been replaced by 3, and -1.23456*10^12 by -3. That is what I called crude!
I have no answer to your second question.

I replaced all floating point numbers by 3 or -3 depending on their signs. As an example, if the number 1.23456*10^12 appeared in your system it would have been replaced by 3, and -1.23456*10^12 by -3. That is what I called crude!
I have no answer to your second question.

I know nothing about MapleSim, but have you looked at
?DAE
in Maple itself?

Since
restart;
factor(expand(3 - 4*cos(z) + cos(2*z))) ;
# returns 2*(cos(z)-1)^2
the inequality holds for all real z.
Thus we get an even more impressive version:
restart:
assume(z>0, beta>0,t::real);
b:= 3*beta*Pi*(3-4*cos(t)+cos(2*t))/(16*z):
-b+sqrt(b^2);

Since
restart;
factor(expand(3 - 4*cos(z) + cos(2*z))) ;
# returns 2*(cos(z)-1)^2
the inequality holds for all real z.
Thus we get an even more impressive version:
restart:
assume(z>0, beta>0,t::real);
b:= 3*beta*Pi*(3-4*cos(t)+cos(2*t))/(16*z):
-b+sqrt(b^2);