The plot clearly shows: one maximum, two minimums and three saddles

Download zadelpunten_procedure_new.mw

Use the  seq  command for this.
Example:

restart;
L:=[0, Pi/6, Pi/4, Pi/3, Pi/2];
[seq([p,sin(p)], p=L)];

restart;
P:=x^4 - x^2 + x - 1:
A:=PolynomialTools:-Split(P, x);
convert(A, radical);
simplify(expand(%));  # A check

1. For  isolve a system should be specified as a set (not a list).
2. You need to specify exactly what you mean by the smallest positive solution, because the solution will not be one number, but an ordered set of numbers (a list). For example, which is less: [2, 4] or [1, 5] ?

In the solution below, the solutions are sorted by the minimum sum of the squares of all coordinates.

restart;
sys:= {a[1] = 2 * a[5], a[1] = a[4], 4 * a[1] + 2 * a[3] = a[6] + 2 * a[7], 2 * a[2] + 2 * a[3] = 2 * a[6], a[2] = a[4] + a[5], 2 * a[2] -a[1] = 2 * a[4]};
indets(sys);
n:=nops(%);
isolve(sys);
assign(%);
solve({seq(a[i]>0, i=1..n)});
[seq(seq([seq(a[i],i=1..n)],_Z2=floor(5*_Z1/2)+1..10),_Z1=1..10)];
sort(%, key = (x->add(x[i]^2,i=1..n)));

I think this contradiction is just a rounding error. Maple uses 10 digits by default. Use precise arithmetic (or make  Digits  bigger):

convert(-24.*y^6*(y-1.)^7*(y+1.)^7*(4.+6.*y^8-22.*y^6+30.*y^4-18.*y^2)/(4.*y^18-28.*y^16+84.*y^14+4.-136.*y^12+116.*y^10-24.*y^8-52.*y^6+56.*y^4-24.*y^2)^2, fraction);
subs(y=1+1/100000, %);
evalf(%);

First, we find the equation of the normal at an arbitrary point  P(x1,x1^2)  on the parabola. This normal intersects the axis  Oy  at a point  Q(0,y0) . Imposing the condition that the distance  PQ  = 1, we first find x1, and then the point  Q (the center of desired circle).

restart;
f:=x->x^2:
g:=x->f(x1)-1/D(f)(x1)*(x-x1):
y0:=g(0):
P,Q:=[x1,f(x1)],[0,y0];
solve((P[1]-Q[1])^2+(P[2]-Q[2])^2=1);
x1:=%[1]:
Q;

# Visualization
A:=plot([f(x),[t,g(t),t=Q[1]..P[1]]], x=-2..2):
B:=plot([P,Q],style=point):
C:=plots:-textplot([[P[],"P"],[Q[],"Q"]], align=right, font=[TIMES,14]):
c:=plottools:-circle(Q,1):
plots:-display(A,B,C,c, scaling=constrained);

Yesterday I answered the same question only asked in another post (this post disappeared for some reason). Here's my answer:
 

restart;
A__left:=1/C__AO-1/C__AA;
f:=1/C__AO;
A__right:=``(f)*(expand(A__left/f));
expand(A__right); # Check

restart;
A__left:=1/C__AO-1/C__AA;
f:=1/C__AO;
A__right:=``(f)*(expand(A__left/f));
expand(A__right); # Check

@jrive  For arbitrary real parameters the last line should be

Pabang := evalc(argument(Pab));

Returns 
                   

From help  "For real arguments x, y, the two-argument function arctan(y, x), computes the principal value of the argument of the complex number 
                            x + I y
, so 
           -Pi < arctan(y, x) and arctan(y, x) <= Pi"

The code below implements this process in Maple. As already pointed out by vv, on each new launch, the numbers  N  may be different.

restart;
randomize():
n:=2012:
S:={seq([i,i], i=1..n)}:
for step from 1 to n-2 do
t:=combinat:-randcomb(S, 2);
S:=(S minus t) union {[n+step,2*(t[1][2]+t[2][2])]};
od:
N:=S[1][2]+S[2][2];
N/1000;

See the thread  https://www.mapleprimes.com/questions/232668-Error-In-Maple-Code

The problem is that Maple does not implement work with piecewise-functions where the number of variables is greater than 1. For example, Maple cannot convert the sum of such functions into some new piecewise-function. Therefore, we have to do something manually. You have 2 variables  lambda[1]  and  lambda2 . Straight lines  (1 - alpha)*lambda[1] - lambda[2]=0  and  -(1 + alpha)*lambda[1] + lambda[2]=0  divide the entire plane  (lambda[1], lambda2)  into 4 regions  R1, R2, R3, R4. In your example,  alpha=0.01 is very small and therefore the regions  R1  and  R3  will be very narrow. Therefore, for better drawing, I took  alpha=0.3 .We define the functions  U1  and  U2  in terms of these regions, and we define the operations on the functions using procedures (for example for `-` sign we use the procedure  `p-`  below) :
 

Download piecewise.mw

restart;
sol:=msolve(x^2=-1, 5);
map2(eval, x, [sol]);

                                    sol := {x = 2}, {x = 3}
                                              [2, 3]

restart;
A:=[seq(a||i,i=1..4)];
B:=[seq(b||i,i=1..4)];
Matrix(4,8, {seq((i,2*i-1)=A[i], i=1..4), seq((i,2*i)=B[i], i=1..4)});

Use the continuous option:

restart;
F:=int( 1/(1+x^n), x=0..1, continuous);
convert(F, hypergeom);