@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.

@aininabdul  Also you missed the multiplication sign after  exp((1/2)*x) . I did not notice this at first.

@Nemo_ 

Try 2 things:

1. Execute  restart  command before the code.
2. Upload here your worksheet or the full text of the code. 

@ThU   in  Maple 2017.3 Windows 32 bit

@asa12 

1. You wrote "is it possible sqrt(3)/6 and 1/2 combine into one term and the rest is one matrix result?". It's impossible.

2. I do not have Maple 12. I did this in Maple 2017.3

Maple by default chooses pretty close colors for graphs, so for greater clarity I suggest setting colors by yourself. scaling  and  size  options are also very useful:

exprs:=[sin(x), cos(x), sqrt(x)*sin(x)];
plot(exprs, x=0..6, color=[red,blue,green], scaling=constrained, legend=exprs, size=[700,400]);

@mrmusic1994  k>0  follows from the physical meaning of the problem, because the resistance force is opposite to the direction of motion.

@tomleslie Thank you for your observation, which allowed me to correctly understand the condition of the problem. At first I understood "right or left" relative to the original coordinate system. In my understanding, this is similar to the movement of a drunken man, who is moving steadily forward, but at the same time he randomly swings to the right or to the left.
I added another option, similar to your one, to my answer.