@vv   Very witty and beautiful solution, vote up! 

It's a pity that the symmetry is broken for any even N.

@vv  Problem point is that each line should intersect with each line. Only in this way we get the maximum number of pieces. I do not know whether there is such a turn-symmetric structure if  n>5

@Mac Dude  Of course this is not the completed code, but only the beginning. OP wants to solve the original partial differential equation by difference method and, therefore, replaces each derivative by the corresponding difference quotient.

@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>]);