@jga  It may be helpful for you that your function can be decomposed into partial fractions:

f:=1/(1+q)^4/(1-q^2)^2*q^4*(q^2+1)*(((q^6+q^5-q^4-2*q^3+q+1)^6/q^13+(q^6+q^5-2*q^3-q^2+q+1)^6/q^25)):
convert(f, parfrac);

@mehdibaghaee  To better understand what you have got, execute at the end  a few commands, for example

A;
A(1, 1);
A(1, 1)[1, 1, 1];
A(1, 1)[2, 2, 5];

@mehdibaghaee   I do not see any errors.

@rlopez  Thank you for reminding about  plots:-shadebetween  command. It appears only in Maple 2015 and I just did not know about this.

@mostafajani  For such functions better use  piecewise  command.

restart; 
f := proc(x) if is(x < 0) then x else x^2 end if end proc;
g := piecewise(x<0, x, x^2);

Examples of use:

plot('f'(x), x=-1..2);
plot(g(x), x=-1..2);
 

@Rouben Rostamian  Strange but somehow I forgot about  this easiest way.

@tomleslie  I just noticed that the equation of the second circle in the initial post was specified with typos, so my answer and the Tom's one  are different.

@vv  Your procedure works well for polynomials with numerical coefficients, but if a polynomial contains parameters, then the procedure fails, for example
GetCoeff(x^2-a*x+5, x);

                                                    0
 

@vv   Perhaps to avoid this, it suffices to require that  f  be homeomorphic mapping?

@Markiyan Hirnyk    Of course, this example works properly only with an affine map. For non-linear mapping every facet of cuboid should be set separately.

restart;
with(plottools): with(plots):
f:=unapply([x^2+3,x*z+3, y-z+3], (x,y,z)); # A nonlinear mapping 
F:=transform(f):
S1:=plot3d([x,0,z], x=0..1, z=0..2.5, style=surface, color=khaki):
S2:=plot3d([x,2,z], x=0..1, z=0..2.5, style=surface, color=khaki):
S3:=plot3d([x,y,0], x=0..1, y=0..2, style=surface, color=khaki):
S4:=plot3d([x,y,2.5], x=0..1, y=0..2, style=surface, color=khaki):
S5:=plot3d([0,y,z], y=0..2, z=0..2.5, style=surface, color=khaki):
S6:=plot3d([1,y,z], y=0..2, z=0..2.5, style=surface, color=khaki):
p:=display(S1,S2,S3,S4,S5,S6):
display(p, F(p), scaling=constrained, axes=normal, orientation=[-15,80], labels=[x,y,z]);

@toandhsp  Replace  %  by  %/sign(lcoeff(%))

@vv   Try this example

PolyOnUnitSphere([ <1,0,0>, <1,2,1>, <0,1,0>, <0,0,1>]);
 

@vv   Very good solution, vote up. But can this technique be useful for a quadrilateral, etc ?

@dantevergil  The shortest ways Acer suggested above. If you do not want to use  remove  or  select  commands, here's another way without packages:

phephoi1 := proc (A, B)
local k, i, C;
k:=0;
for i from 1 to nops(A) do
   if A[i] in B then k:=k+1; C[k]:=A[i]; fi;
end do;
convert(C, list);
end proc:

Example of use:

A := [1, 2, 3, 5, 6, 7];
B := [2, 4, 7, 11, 8];
phephoi1(A, B);

@acer   Oh sure. My goal was to keep the idea of initial code, but to make it work.