See the  Large Operators  palette :

In Maple we can use the  extrema  command for symbolic calculation of the  maximum  or the  minimum . The  extrema  command employs the method of Lagrange multipliers. Of course, this method is not universal and can fail if the extremum is reached at the points where differentiability is violated.

We see that the maximum  sqrt(2)/4  is reached at two points. At the same time, we also learned everything about the minimum.

Download extrema.mw

1. For grouping the terms of an equation in Maple do not use square or curly brackets (only parentheses).

2. There is no need to use the  fsolve  command first for plotting, you can immediately plot with the  plots:-implicitplot  command. We see 2 values of  M  for each value of  sigmai from the range  0.1..0.4 :

Download eq_new.mw

Everything you need to write the procedure can be easily found in these two articles:
https://en.wikipedia.org/wiki/Conic_section#:~:text=There%20are%20three%20types%20of,considered%20as%20a%20fourth%20type.&text=If%20the%20cutting%20plane%20is,and%20is%20called%20a%20parabola.
and
https://en.wikipedia.org/wiki/Degenerate_conic

We can think of a conformal mapping as a mapping from R^2 to R^2 :

>  restart; assume(y, real); assume(x, real); f := z -> I + z*exp(1/4*I*Pi); w := f(x + y*I); u := Re(w); v := Im(w); with(plots): with(plottools): R:=display(polygon([[0,0],[1,0],[0,1]],color="LightBlue"),seq(plot(i,x=0..1-i,color=black),i=0..1,0.1),plot([seq([i,t,t=0..1-i],i=0..1,0.1)],color=black)): F:=transform(unapply([u,v],x,y)): `f(R)`:=display(F(R)): display(<R |`f(R)`>, scaling=constrained, size=[300,300]); display(R,`f(R)`); >

Download conf.mw

Download inequal_question_new.mw

restart;
eq:=x^2+2*x-1=0;
x=~[solve(eq)];

Let  a  be the number of points on one dice, and   b  be the number of points on another one.
 

Download LE.mw

You can use an inert form for subtraction to keep the expression unchanged. For any calculations and transformations use the  value  command:

restart;
expr:=arccos(p%-a);
value(expr);
value(eval(expr,[p=0,a=1]));

_B1  is a binary variable that is 0 or 1 :

restart;
ode:=diff(y(x),x)-y(x)/x+csc(y(x)/x)=0;
sol:=dsolve([ode,y(1)=0]);
about(_B1);
sol1:=eval(sol,_B1=0);
sol2:=eval(sol,_B1=1);
odetest(sol1,ode);
odetest(sol2,ode);

Here is a solution in Maple:
 

Download nlsystem.mw

restart;
plots:-animate(plot,[cos(2*theta),theta=0..a, coords=polar], a=0..2*Pi);

We can significantly increase efficiency if we use symmetries: from one solution we can get  6 * 6 = 36 solutions by permutations of the first three and last three numbers. Instead of loops, it’s slightly more efficient and more compact to use nested seq . So we get the unique solution:

Download isys.mw

Your surface is a torus:

restart;
plot3d(eval([x,sqrt(x^2+y^2)*cos(s),sqrt(x^2+y^2)*sin(s)],[x=cos(t),y=2+sin(t)]), t=0..2*Pi, s=0..2*Pi, scaling=constrained, labels=[x,y,z]);

Download shade.mw

Here is another shorter way:
 

Download shade1.mw

Here is a more automatic plotting for the intersection:

restart;
r1:=2:  r2:=2*(1-cos(theta)):
plot(min(r1,r2), theta=0..2*Pi, color=green, coords=polar, filled, scaling=constrained);

Edit.