f:= x-> (sign(4*x^2-4)*abs(4*x^2-4)^(1/5))^4;
evalf(f(0));

Use  plots:-display  command for this:

plots:-display(rect[i], i=1..21);

Edit: of course above should be
plots:-display(seq(rect[i], i=1..21));

Thanks to Carl for the amendment.
 

You can treat this as a curve in 3d.

Example:

restart;
eq1:=sin(z):
eq2:=cos(z):
plots:-spacecurve([eq1,eq2,z], z=0..4*Pi, color=red, thickness=3);

Everything works with assuming:

simplify(A) assuming real;
simplify(B) assuming real;
simplify(A) assuming complex;
simplify(B) assuming complex;
combine(A);
combine(B);
combine(A*B);

                                                          1
                                                          1
                                                          1
                                                          1
                                                          1
                                                          1
                                                          1

Or

simplify(numer(A))/denom(A);

and so on ...

Several other ways:

is(A=1);
is(B=1);
is(A*B=1);
combine(convert(A, tan));
combine(convert(B, tan));
combine(convert(A*B, tan));

Edit.

I swapped  seq  and  plottools:-arrow  and added a few required parameters for  plottools:-arrow  command.

plot_arrows_sequence_new.mw

Addition.  [seq] means that the result of applying  seq  command will be surrounded by square brackets, i.e. we get a list.

Example:

seq(k^2, k=1..10);
[seq](k^2, k=1..10);
        

You can do this using  LinearAlgebra:-GenerateMatrix  command as below:

Download MatrixForm1.mw

For matrices and vectors use  LinearAlgebra:-Equal  command:

for i from 2 while not LinearAlgebra:-Equal(M[i, 1], M[1, 1]) and i < 25 do
print(i, M[i, 1]); 
end do;

2 errors were fixed:
1. gamma is a protected constant in Maple. So I wrote  local gamma;
2. deq1  and   deq2  should be  instead of  deg1  and  deg2

Addition - the plottings:

odeplot(R, [[eta,f(eta)],[eta,theta(eta)],[eta,phi(eta)]], eta=0..8, legend=[f(eta),theta(eta),phi(eta)]);
      

Download dsolve.mw

Here I try to explain why  plots:-implicitplot  command (without special methods like acer' ones) does not work. If I understand correctly, it works simply by calculating the values ​​of the function on the grid without using  fsolve  and the like. If at two points  (x1,y1)  and  (x2,y2)  the continuous function  F(x,y)  takes values ​​of different signs, then at some point of the segment with the ends   (x1,y1)  and  (x2,y2) it is  0 . In the first very simple example below  (x-sqrt(2))^2=0,  implicitplot  does not plot anything, because there is no change of signs.

Next, we find the domain of the original function  f(b,p) , using the fact that under the roots there are homogeneous functions. As a result, we find that this domain of definition is located between two very close straight lines  p=0.2466777778*b  and  p=0.25*b .

Since the curve  F(b,p)=0  itself is very close to the lower boundary of this domain, if we plot in large enough ranges (for example for  b=0..5, p=0..5) Maple simply skips the area of the sign change, because the grid sizes are much larger. Only if we consider a small rectangle (I took the range  b = 7.95..8.05, p = 1.95..2, the plotting is happening (see the bottom 2 graphics) :

Download Plottings.mw

Here is another method based on applying  extrema  command and simplifying the result by  rationalize  and  simplify  commands:

restart;
f:=x->(cos(x)+sqrt(3)*sin(x))/(cos(x)+sin(x)+2);
convert(extrema(f(x), {}, x), list);
Min, Max:=simplify(rationalize(%))[];
evalf([Min,Max]);

Output:

You assigned  o8:=0  and  o8  stands in the denominators.

In Maple 2018 everything is the same.
Do not use  simplify  command, but just differentiate:

with(Physics):
Setup(noncommutativeprefix={P,Q});
diff(Q(t)^2*P(t)*Q(t) + Q(t)*P(t)*Q(t)^2, t);

Output:

Replace  e^(c*v)  by  exp(c*v) . e  is just a symbol in Maple.

This works:

restart:
with(Statistics): 
X := Vector([1, 2, 3, 4, 5, 6], datatype = float):
Y := Vector([2.2, 3, 4.8, 10.2, 24.5, 75.0], datatype = float): 
NonlinearFit(a+b*v+exp(c*v), X, Y, v);

sqrt(sqrt(p^2+1)-1)*sqrt(sqrt(p^2+1)+1);
simplify(combine(%)) assuming positive;

Here is another rather effective method of simplifying expressions containing radicals. Since this expression itself (we denote it by  A ) is obviously positive for any real  p , it is equal to the square root of its square   A=sqrt(A^2) . Therefore, we first square it, then simplify and extract the square root:

A:=sqrt(sqrt(p^2+1)-1)*sqrt(sqrt(p^2+1)+1);
sqrt(simplify(A^2))  assuming p>0;
                         

  Edit.  

If you define a matrix or an array, you should specify ranges for the indexes and just specify non-zero elements. The rest of the elements are automatically considered zeros. 

Your example:

H:=Array((1..3)$4, {(1,1,1,1)= value1, (1,2,2,1) = value2, ...});

See help on  Matrix  and  Array  commands for details.