@Jjjones98 

restart;

y:= x-> ((1 + a*x + b*x^2)/(1 + c*x + d*x^2))*(ln(sinh(x)^2 + cosh(x)^2)):

Sys:={seq(coeff(series(y(x), x=0, 7), x^n)=0, n=3..6)}:

Sol:=solve(Sys);

series(eval(y(x),Sol), x=0, 7);

P:=convert(%, polynom);

eval(y(x), Sol)=P;

solve(%, ln(sinh(x)^2 + cosh(x)^2));

The questions remain.

Why it works:
patmatch(x^2*y^(-2), x^n::realcons*y^m::realcons, 's');  
s;
                              true
                        [m = -2, n = 2]

But it does not work:
patmatch(x^2/y^2, x^n::realcons/y^m::realcons, 's');
s;
                             false
                               s


 

@vv  I completely agree. I never use  parametric  option also because Maple is very weak in working with parameters. Tasks with parameters I always solve with Mathematica by Reduce command.

@animeplot  I do not have Maple 13 and so I can not check it. In Maple 2015 - 2017 everything works correctly. 

Try the following code. I think that it will work:

restart;
L1:=s->[[0,0],[5,0],[2.5,s],[0,0]]:
L2:=s->[[2,0],[7,0],[4.5,s],[2,0]]:
L3:=s->[[4,0],[9,0],[6.5,s],[4,0]]:
A:=s->plot([L1(4-round(s)),L2(4-round(s)),L3(4-round(s))], color=red, scaling=constrained, axes=none):
N:=15:  
plots:-animate(A, [n], n=0..N, frames=75, background=A(0));

Edit.

@animeplot  Yes we can do it. Using  plots:-animate  command is usually preferable, the code is easier. To make the animation look discrete, I used  round(n)  instead  n . Here is your code using plots:-animate  command:

restart;
L1:=s->[[0,0],[5,0],[2.5,s],[0,0]]:
L2:=s->[[2,0],[7,0],[4.5,s],[2,0]]:
L3:=s->[[4,0],[9,0],[6.5,s],[4,0]]:
A:=s->plot([L1(s),L2(s),L3(s)], color=red, scaling=constrained, axes=none):
N:=15:  
plots:-animate(plots:-display,[A(4-round(n))], n=0..N, frames=75, background=A(4));

@ANANDMUNAGALA   Thank you. Write in more detail for which bodies you want to see similar cross-sections.

@Les  It seems that you misunderstood me. I meant that one of the possible ways to solve your problem will be the use of recursion. A recursive procedure is a procedure that in its code refers to itself. For example, the factorial of a positive integer can be found using a simple recursive procedure:

F:=proc(n)
if n=1 then return 1 fi;
n*F(n-1);
end proc:
F(10);
                              3628800

Of course, for your problem the solution will be much more difficult.

@nm  Thank you. I did not consider it proper to send my reply as an answer (because the main idea is not mine), but at your request I did it.

Addition. Here are the answers to your 2 questions:

1. 
expr:=(x^2*y+y)*y^3*exp(-x-y+1)*3^(-x-y)*sqrt(x^2*y-y):
r:=Splitting(expr):
xPart:=op(r)[1];           
yPart:=op(r)[2];
simplify(expand(xPart*yPart-expr), symbolic);

2. You wrote about an arror for  expr:= -Pi*(1-cos(x)^2)*sin(y)^3/exp(1);

The reason for the error is that Maple for some reason treats the constant  Pi  as a name. Here is a small adjustment procedure that removes this problem:

Splitting1:=proc(Expr)
local x, y, Expr1, Expr2;
x,y:=(indets(Expr, name) minus {constants})[];
Expr1, Expr2 := selectremove(z -> has(z,x) and has(z,y), simplify(factor(Expr), symbolic));
``(simplify(select(has, Expr2*expand(Expr1), x)))*
``(simplify(remove(has, Expr2*expand(Expr1), x)));
end proc:

expr:= -Pi*(1-cos(x)^2)*sin(y)^3/exp(1);
r:=Splitting1(expr);  # OK

​​​​​​​Edit.
 

By analogy with the formula  n!=n*(n-1)! , your problem can easily be solved by means of a recursive procedure for an arbitrary  n . The resulting graphs can be conveniently located in the form of an array consisting of  n  individual trees.

@John Fredsted  factor  command is necessary. Consider the following example

Expr:=x*y-sqrt(2)*x+sqrt(2)*y-2;
 

It is convenient to use  `if` construction here. print command can be omitted:

a := 1:
b := 3:
`if`(a>b, wrong, correct);
                                                  correct

  @sajad
OK, I just did not notice the exponent of degree 2. Here is the corrected solution:

System_new.mw

Explanation. The idea of the solution is the same as in my previous version, that is, the elimination of the integral. But this has to be done 2 times, because we have two different integrals and we get the differential equation of order 2. Therefore, we have to look for 2 initial conditions.

@acer  
restart;
simplify(solve({x=r*cos(theta), y=r*sin(theta), r>0}, {r,theta}, explicit))  assuming real;

@Ramakrishnan  I can not understand the reasons for the errors, because I do not have enough information. What is  BP  and what are the elements of the matrix  M1. Here's what I see when I try to run your code:

@Ramakrishnan 

The command  plot([[x1,y1],[x2,y2]], x=a..b)  just builds the segment that connects the points  [x1,y1]  and  [x2,y2] , regardless of the view  x=a..b

But if we want to continue this segment to the whole range a..b , then we need the equation of the corresponding straight line and  g  procedure provides us with this equation. x1, y1, x2, y2 are local variables in this procedure and their names do not matter. For example, g:=(a,b,c,d)->b+(d-b)*(x-a)/(c-a)  will be the same.