I do not understand what the problem is. Maple handles this task easily:

Download rec.mw

Unfortunately, Optimization:-Maximize command in this example returns an erroneous result (I use Maple 2018.2), since it returns only a local extremum. For the correct solution, it is useful to plot graphs and use  initialpoint  option:
 

Edit.

Download Max1.mw

As a workaround use the  eval  command:

restart;
f(x):=a+b;
a:=5;
b:=2;
f(x):=eval(f(x));

Edit. In addition to the acer's advice, you can also use indexed names:

restart;
f[x]:=a+b;
a:=5;
b:=2;
f[x]:=f[x];

gamma  is a protected constant in Maple, and  S  cannot be differentiated by a constant. Execute it first

local gamma;

We have

f = 2*x^5-x^3*y+2*x^2*y^2-x*y^3+2*y^5 = 2*(x^5+y^5) - x*y*(x^2+y^2-2*x*y) = 2*(x^5+y^5) - x*y*(x-y)^2

If  x<0  and  y<0  then both summands are negative.

A simple procedure called  `lim` allows you to get the result in the desired form. I did not find in Maple the possibility of directly exporting this result to MathType, but you can first convert it to latex, and then find the conversion of latex to MathType on the Internet:

Download conversions.mw

You can easily verify this identity for any  n :

restart;
F:=unapply(rsolve({F(n+1) = F(n)+F(n-1),F(1)=1,F(2)=1},F(n)),n);
combine(expand(F(n+1)*F(n+2) - F(n)*F(n+3)));

Edit. The proof by mathematical induction:

restart;
simplify(subs(n=k+1,F(n+1)*F(n+2) - F(n)*F(n+3)), {F(k+1)*F(k+2) - F(k)*F(k+3)=a,F(k+2)=F(k)+F(k+1),F(k+3)=F(k+1)+F(k+2),F(k+4)=F(k+2)+F(k+3)});
eval(%, a=(-1)^k);
simplify(%-(-1)^(k+1));

f:=(i,j)->``((x-GaußKnoten[j])/(GaußKnoten[i]-GaußKnoten[j]));
n:=3;
`*`(seq(seq(`if`(i<>j,f(i,j),NULL), j=1..n), i=1..n));

    
Edit. In order to explicitly multiply all this, use the  expand  command (code continuation):

expand(%);

Using the  plots:-odeplot  command, we can easily not only plot the curve, but also animate it:

plots:-odeplot(F, [x(t), y(t), z(t)], t = 0 .. 0.7, color = red, thickness = 2, axes = normal, scaling = constrained)

plots:-odeplot(F, [x(t), y(t), z(t)], t = 0 .. 0.7, color = red, thickness = 2, axes = normal, scaling = constrained, frames = 15);

Maybe the following more detailed form of writing will be clearer:

diff(w(y(t), t), t);
convert(%, diff);

This is easy to do if you use the plots:-inequal  command and Cartesian coordinates:

with(plots):
P1 := plot([-sin(t), t, t = 0 .. 2*Pi], coords = polar, color = red):
P2 := plot([cos(t), t, t = 0 .. 2*Pi], coords = polar, color = blue):
Shade := inequal({(x-1/2)^2+y^2<1/4,x^2+(y+1/2)^2<1/4}, x=-0.5..1, y=-1..0.5, color=yellow, nolines):
display(P1, P2, Shade, scaling = constrained);

Below is another method based on approximating a region on the plane by a polygon: 

with(plots): with(plottools):
P1 := plot(-sin(t), t = 0 .. 2*Pi, coords = polar, color = red):
P2 := plot(cos(t), t = 0 .. 2*Pi, coords = polar, color = blue):
P:=polygon([seq([-sin(t)*cos(t),-sin(t)*sin(t)],t=-evalf(Pi/4)..0,0.1),seq([cos(t)*cos(t),cos(t)*sin(t)],t=-evalf(Pi/2)..-evalf(Pi/4),0.1)],color=yellow,style=surface):
display(P1, P2, P, scaling = constrained);

This method is automated in my post  https://www.mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-  and is applicable to any flat region bounded by a piecewise smooth non-selfintersecting curve (the curve can be specified in Cartesian or polar coordinates or parametrically).

Edit.

As alternatives to  subs , you can also use  eval  or  alias :

restart;
s1 := RootOf(D1*D2*D6*_Z^2+(-D1*D4*D6-D2*D6*S-D4)*_Z+D4*D6*S):

s2 := -(D1*RootOf(D1*D2*D6*_Z^2+(-D1*D4*D6-D2*D6*S-D4)*_Z+D4*D6*S)*D6-S*D6+RootOf(D1*D2*D6*_Z^2+(-D1*D4*D6-D2*D6*S-D4)*_Z+D4*D6*S))/D2/D6/RootOf(D1*D2*D6*_Z^2+(-D1*D4*D6-D2*D6*S-D4)*_Z+D4*D6*S):

eval(s2, s1=s);
alias(s=s1):
s2;

If  subs  and  eval  are fairly close commands, then  alias  is essentially different from them. If in the first two cases  s  is just a symbol, then as for  alias  Maple remembers what  s  means. This may be more convenient in subsequent calculations. See help on these commands.

restart;
lambda:=x+I*y:
evalc(abs(1+1/3*lambda+1/18*lambda^2-1/324*lambda^3+1/1944*lambda^4)-1);
plots:-implicitplot(%, x=-15..15, y=-15..15, gridrefine=5, scaling=constrained, size=[550,550]);

Edit. It would be interesting to show that the curves on the right (bottom and top) are exact circles (this is my assumption) and learn their equations.

The following is working:

Download plot.mw

Download Maple_Inverse_Matrix_test_new.mw