Boundary-layer ODE BVPs like this related to nano-fluid dynamics are one of the most common sources of Questions here on MaplePrimes (perhaps even the single largest source), and I've worked through several such papers here. In particular, see my Post "Numerically solving BVPs that have many parameters" . Here is how to read this paper with an eye towards computing the solutions in Maple. The parts of the paper that you actually need to use are extremely short.

1. Everything in the first two sections (10 pages) of the paper can be ignored, except to the extent that after the solution procedure is coded, we'll need to look in those sections for the ranges of numeric values of the parameters to use.

2. The complete BVP such as it's needed for computation is given in Eqs. 21 - 23 on page 11. That's all there is! Here, the authors have renamed the dependent variables of the system as y(1), ..., y(6), yy__1, yy__2. That's a complete load of crap! Just stick with the original names f(eta), ..., diff(f(eta), eta$4), theta(eta), ..., diff(theta(eta), eta$2), which are the names commonly used in nano-fluid dynamics (as I'm sure you're aware, this being not the first time that you've posted to MaplePrimes about such a BVP). If the authors have done this renaming because Matlab syntax requires it, then shame, shame on Matlab!, because this notation is hideous. Maple handles this conversion to a first-order system in the background.

All parameters that appear symbolically in Eqs. 21 - 23 must remain symbolic in your coding of the BVP at this point. That is, we won't assign them numeric values until dsolve is called.

3. One of the boundary conditions is stated "as eta -> infinity". Maple's BVP solver doesn't allow for that. We simply pick a value (such as 1) to use instead of infinity. But don't directly use 1! That would be a very bad programming decision, likely to lead to confusion, error, and/or inflexibility. Instead, give this "fake infinity" a name (such as inf), and use that exclusively for the rest of the code. Thus, that 6th boundary condition should be theta(inf) = 0. It'll be worthwhile to look through the paper more thoroughly to see how the authors handled this. 

4. Everything up until this point is extremely straightforward, and the details have been shown hundreds of times (at least) on MaplePrimes.

5. The major hurdle with these boundary-layer BVPs (and the cause of a great number of those MaplePrimes Questions) is a failure to converge for some values of the parameters due to non-uniqueness of the solution. (In my limited experience, this is more likely to happen when a "fake infinity" is used.) A method to handle this is discussed in the first paragraph of section 4. The syntax for doing this in Maple uses the approxsoln option to dsolve.

It will be interesting to see if Maple can perform as well as, or better than, Matlab for this. From a computational viewpoint, this is the only interesting part of the problem, because, as I said, everything else is completely straightforward. But this final part is quite interesting.

The entirety of this process can be done without the slightest bit of knowledge of fluid dynamics; indeed, I have no such knowledge myself. It's a purely a Numerical Analysis problem.

@Wilks Yes, what you just said about restart is correct, and thus you can't directly use Kitonum's Answer above. I think that he was assuming that Y1 and its derivative were undefined functions.

Assuming that your `Y1"`(x) is a correctly defined function (or algebraic expression) for which Maple can find a symbolic antiderivative, your C1 can be expressed simply as 

-intat(`Y1"`(x), x= 0)

If you prefer, you could also write explicitly C1 = -intat(`Y1"`(x), x= 0). Regardless of the complexity, there'd never be a need to use solve for this; however, if for some reason (pedagogy, perhaps) you'd prefer to use solve, you could do

solve({intat(`Y1"`(x), x= 0) + C1}, C1)

If you want a procedure to do this for arbitrary initial conditions, you could use this:

IC:= proc(`y'`::appliable, ic::function(algebraic)=algebraic, C::name)
local x;
    `tools/genglobal`(`if`(nargs=3, C, ':-_C')) = 
        rhs(ic) - intat(`y'`(x), x= op(lhs(ic)))
end proc:

And you'd use it like this (assuming that `Y1"` has been defined, Y1 has not been defined,  x0 and y0 are any appropriate algebraic expressions (including simply numbers), and C1 is an unassigned name):

IC(`Y1"`, Y1(x0)=y0, C1);

where the 3rd argument, C1, is optional, and if not supplied the procedure will choose an unused name beginning with _C.

The command forget can be used to dump the cache tables and remember tables. See its help page ?forget. However, in my opinion (based on many years experience intensively time testing Maple code), more meaningful results can be obtained by calling the procedure multiple times with random input. Generate a random list of sets of polynomials (of the type that you're interested in), and map the procedure over that list.

Use CodeTools:-Usage, and note both the real time and cputime. 

Note that your desired unexpanded form f^6/6 has a constant term of 2048/3. This constant term doesn't appear in Maple's int result, so that result can't be factored to f^6/6.

The immediate first thing that I noticed was the same thing discussed by acer. But there are several other issues, and I decided to show you how to do a more thorough analysis of your code.

These are some steps to take that will quickly show you which line of code is taking so long. But first I did some preliminary steps to avoid some obviously unnecessary computations:

• I replaced all int(...) with evalf(Int(...)). Note the use of capital I!
• I replaced all sum with add.

Next, I dealt with your calls to pprint. I assume that this is a procedure that you use to print supplementary information perhaps as an aid in debugging. The problem with that is that even if the procedure does nothing at all (indeed, even if the procedure doesn't even exist), it still takes a significant amount of time simply to evaluate its arguments. There is a built-in command made just for this situation: userinfo. I replaced all pprint(...) with userinfo(1, Erweiterung, ...). This call will evaluate and print the arguments if and only if infolevel[Erweiterung] >= 1 where infolevel is a global table that's already predefined explicitly for use with userinfo. (It's traditional to use a procedure's name as the 2nd argument to userinfo, but actually any name can be used. The name must match the infolevel entry.)

Next, I used the command CodeTools:-Profiling:-Profile to "profile" Erweiterung. "Profiling" means to create a line-by-line analysis showing the number of times each line is executed and how much time and memory it used. This is cumulative over all calls to the procedure in the current session; so, if you want to analyze a single call to the procedure, start with restart. Thus:

#The procedure code and (very importantly!) the restart are above. 
#I just didn't want to clutter up this Answer with that code.

Digits:= 15: #experimental value

(*
#The following line of code is commented out because I *don't* 
#want the userinfo calls to take any time.
infolevel[Erweiterung]:= 1:
*)

CodeTools:-Profiling:-Profile(Erweiterung):

#This was your troublesome call:
Erweiterung(-1,1,2*x^2+2,2/sqrt(1-x^2),[-.8660254037, 0, .8660254037], 4);

DieKnoten, -0.999999999950535, -0.866025403700000, 
  -0.499999999951372, 0., 0.499999999951372, 0.866025403700000, 
  0.999999999950535

DieGewichte, TABLE([1 = 0.523598776221123, 2 = 1.04719755075381, 
  3 = 1.04719755105594, 4 = 1.04719755111737, 
  5 = 1.04719755105594, 6 = 1.04719755075381, 
  7 = 0.523598776221123])

                        18.8495559215366

To view the "profile", use CodeTools:-Profiling:-PrintProfiles:

CodeTools:-Profiling:-PrintProfiles(Erweiterung);
Erweiterung
Erweiterung := proc(Unten, Oben, f, g, Liste, n)
local Basenwechsel, p, i, c, d, e, Hn, Koeffizienten, a, s, j, M, V, S, K, nNeu, Em,
Hnm, KnotenHnm, KoeffizientenHnm, h0, b, gxi, Gewichte, Delta, Ergebnis, 
  Endergebnis, Koeffizient, Rest, Basenwechsel1, Basenwechsel2, k;
     |Calls Seconds  Words|
PROC |    1   1.438 17286788|
   1 |    1   0.000     37| Basenwechsel1 := proc( Dividend, m ) :: list;
                                local Basenwechsel;
                                Koeffizient := quo(Dividend,p[m],x);
                                Rest := rem(Dividend,p[m],x);
                                if m = 0 then
                                    Basenwechsel := [Koeffizient*evalf(Int(g*p
                                      [m]^2,x = Unten .. Oben))];
                                else
                                    Basenwechsel := [Koeffizient*evalf(Int(g*p
                                      [m]^2,x = Unten .. Oben)), op(procname(
                                      Rest,m-1))];
                                end if;
                            end proc;
   2 |    1   0.000     31| Basenwechsel2 := proc( Dividend, m ) :: list;
                                local Basenwechsel;
                                Koeffizient := quo(Dividend,p[m],x);
                                Rest := rem(Dividend,p[m],x);
                                if m = 0 then
                                    Basenwechsel := [Koeffizient];
                                else
                                    Basenwechsel := [Koeffizient, op(procname(
                                      Rest,m-1))];
                                end if;
                            end proc;
   3 |    1   0.000    273| p[-1] := 0;
   4 |    1   0.000     10| p[0] := 1;
   5 |    1   0.000      6| for i to 2*numelems(Liste)+2*n do
   6 |   14   1.266 15477600|     p[i] := x^i-add(evalf(Int(x^i*p[j]*g,x = Unten
                                  .. Oben))*p[j]/evalf(Int(p[j]^2*g,x = Unten
                                  .. Oben)),j = 0 .. i-1);
   7 |   14   0.000      0|     userinfo(1,Erweiterung,p[i]);
   8 |   14   0.000    778|     c[i-1] := coeff(p[i],x,i)/coeff(p[i-1],x,i-1);
   9 |   14   0.000   1149|     d[i-1] := (coeff(p[i],x,i-1)-c[i-1]*coeff(p[i-\
                                  1],x,i-2))/coeff(p[i-1],x,i-1);
  10 |   14   0.000     42|     if i <> 1 then
  11 |   13   0.000   2504|         e[i-1] := coeff(p[i]-(x*c[i-1]+d[i-1])*p[i
                                      -1],x,i-2)/coeff(p[i-2],x,i-2)
                                else
  12 |    1   0.000    275|         e[i-1] := 0
                                end if
                            end do;
  13 |    1   0.000      0| userinfo(1,Erweiterung,Liste[1],numelems(Liste));
  14 |    1   0.000     52| Hn := mul(x-Liste[i],i = 1 .. numelems(Liste));
  15 |    1   0.000      0| userinfo(1,Erweiterung,Hn);
  16 |    1   0.015  34905| Koeffizienten := PolynomialTools:-
                              FromCoefficientList(Basenwechsel1(Hn,numelems(
                              Liste)+1),x,termorder = reverse);
  17 |    1   0.000      0| userinfo(1,Erweiterung,Koeffizienten,HIER);
  18 |    1   0.000      0| userinfo(1,Erweiterung,c,d,e);
  19 |    1   0.000    549| a[0][0] := 1;
  20 |    1   0.000    286| a[1][0] := x;
  21 |    1   0.000     68| a[1][1] := -e[1]*c[0]/c[1]+(d[0]-d[1]*c[0]/c[1])*x
                              +c[0]/c[1]*x^2;
  22 |    1   0.000      2| for s from 2 to numelems(Liste)+n do
  23 |    6   0.000   1746|     a[s][0] := x^s;
  24 |    6   0.000    620|     a[s][1] := -e[s]*c[0]/c[s]*x^(s-1)+(d[0]-d[s]*
                                  c[0]/c[s])*x^s+c[0]/c[s]*x^(s+1);
  25 |    6   0.000    194|     pprint(coeff(a[s][1],x,s),as1s)
                            end do;
  26 |    1   0.000      2| for s from 2 to numelems(Liste)+n do
  27 |    6   0.000      0|     for j from 2 to s do
  28 |   21   0.000      0|         userinfo(1,Erweiterung,c[j-1]*add(coeff(a[
                                      s][j-1],x,k-1)/c[k-1]*x^k,k = abs(s-j)+2
                                      .. s+j));
  29 |   21   0.000      0|         userinfo(1,Erweiterung,add((d[j-1]-c[j-1]*
                                      d[k]/c[k])*coeff(a[s][j-1],x,k)*x^k,k =
                                      abs(s-j)+1 .. s+j-1));
  30 |   21   0.000      0|         userinfo(1,Erweiterung,c[j-1]*add(e[k+1]*
                                      coeff(a[s][j-1],x,k+1)/c[k+1]*x^k,k = 
                                      abs(s-j) .. s+j-2));
  31 |   21   0.000      0|         userinfo(1,Erweiterung,e[j-1]*add(coeff(a[
                                      s][j-2],x,k)*x^k,k = s-j+2 .. s+j-2));
  32 |   21   0.000  31718|         a[s][j] := c[j-1]*add(coeff(a[s][j-1],x,k-\
                                      1)/c[k-1]*x^k,k = abs(s-j)+2 .. s+j)+add
                                      ((d[j-1]-c[j-1]*d[k]/c[k])*coeff(a[s][j-\
                                      1],x,k)*x^k,k = abs(s-j)+1 .. s+j-1)-c[j
                                      -1]*add(e[k+1]*coeff(a[s][j-1],x,k+1)/c[
                                      k+1]*x^k,k = abs(s-j) .. s+j-2)+e[j-1]*
                                      add(coeff(a[s][j-2],x,k)*x^k,k = abs(s-j
                                      )+2 .. s+j-2)
                                end do
                            end do;
  33 |    1   0.000      2| for s from 0 to numelems(Liste)+n do
  34 |    8   0.000      0|     for j from 0 to s do
  35 |   36   0.000      0|         userinfo(1,Erweiterung,a[s][j],Polynom[s][
                                      j])
                                end do
                            end do;
  36 |    1   0.000   6724| M := Matrix(n,n);
  37 |    1   0.000   3919| V := Vector(n);
  38 |    1   0.000      0| for s from 0 to n-1 do
  39 |    4   0.000      0|     for j from 0 to s do
  40 |   10   0.000   1674|         M(s+1,j+1) := add(coeff(a[s][j],x,k)*coeff
                                      (Koeffizienten,x,k),k = 0 .. n);
  41 |   10   0.000     30|         if s <> j then
  42 |    6   0.000     54|             M(j+1,s+1) := M(s+1,j+1)
                                    end if;
  43 |   10   0.000      0|         userinfo(1,Erweiterung,M,1)
                                end do;
  44 |    4   0.000      0|     userinfo(1,Erweiterung,testb1);
  45 |    4   0.000      0|     userinfo(1,Erweiterung,coeff(a[n][s],x,2));
  46 |    4   0.000      0|     userinfo(1,Erweiterung,coeff(Koeffizienten,x,2
                                  ));
  47 |    4   0.000      0|     userinfo(1,Erweiterung,testb2);
  48 |    4   0.000      0|     userinfo(1,Erweiterung,Koeffizienten);
  49 |    4   0.000    955|     M(s+1,n+1) := add(coeff(a[n][s],x,k)*coeff(
                                  Koeffizienten,x,k),k = 0 .. n);
  50 |    4   0.000      0|     userinfo(1,Erweiterung,M,V)
                            end do;
  51 |    1   0.000      0| userinfo(1,Erweiterung,M,V);
  52 |    1   0.016 193237| S := LinearAlgebra:-LinearSolve(M,V);
  53 |    1   0.000     48| K := evalindets(S,name,() -> 2);
  54 |    1   0.000      0| userinfo(1,Erweiterung,K,LinSolve);
  55 |    1   0.000    260| Em := add(p[i]*K[i+1],i = 0 .. n);
  56 |    1   0.000     14| Hnm := Hn*Em;
  57 |    1   0.016 161926| KnotenHnm := fsolve(Hnm,complex);
  58 |    1   0.000      0| userinfo(1,Erweiterung,Hn,alt,Em,neu,Hnm);
  59 |    1   0.000      0| userinfo(1,Erweiterung,Testergebnis,nNeu);
  60 |    1   0.000      0| userinfo(1,Erweiterung,fsolve(Hnm),fsolve(nNeu));
  61 |    1   0.000  11230| KoeffizientenHnm := ListTools:-Reverse(
                              Basenwechsel2(Hnm,numelems(Liste)+n));
  62 |    1   0.000      0| userinfo(1,Erweiterung,KoeffizientenHnm);
  63 |    1   0.000     12| h0 := evalf(Int(g,x = Unten .. Oben));
  64 |    1   0.000      0| userinfo(1,Erweiterung,h0,HO);
  65 |    1   0.000    277| b[n+numelems(Liste)+2] := 0;
  66 |    1   0.000     14| b[n+numelems(Liste)+1] := 0;
  67 |    1   0.000      4| for i to nops([KnotenHnm]) do
  68 |    7   0.000     14|     for j from numelems(Liste)+n by -1 to 1 do
  69 |   49   0.000      0|         userinfo(1,Erweiterung,test21,KnotenHnm,
                                      Hnm);
  70 |   49   0.000   2108|         b[j] := KoeffizientenHnm[j+1]+(c[j]*
                                      KnotenHnm[i]+d[j])*b[j+1]+e[j+1]*b[j+2];
  71 |   49   0.000      0|         userinfo(1,Erweiterung,b[j])
                                end do;
  72 |    7   0.000      0|     userinfo(1,Erweiterung,test23);
  73 |    7   0.000  21442|     gxi := quo(Hnm,x-KnotenHnm[i],x);
  74 |    7   0.000      0|     userinfo(1,Erweiterung,gxi);
  75 |    7   0.000   1230|     Gewichte[i] := c[0]*b[1]*h0/eval(gxi,x = 
                                  KnotenHnm[i]);
  76 |    7   0.000      0|     userinfo(1,Erweiterung,b);
  77 |    7   0.000    361|     Delta[i] := c[0]*b[1]
                            end do;
  78 |    1   0.125 1293248| print(DieKnoten,KnotenHnm);
  79 |    1   0.000  34818| print(DieGewichte,Gewichte);
  80 |    1   0.000    340| Ergebnis := add(eval(f,x = KnotenHnm[k])*Gewichte[
                              k],k = 1 .. nops([KnotenHnm]));
  81 |    1   0.000      0| Endergebnis := evalf(Ergebnis)
end proc

And so we see that the line of code labelled 6 in the left column is where the vast majority of the the time is being used. Line 78 is also a problem, and likely should be replaced with userinfo. Interpretting the memory numbers is more of an expert-level topic and you shouldn't concern yourself with that at the moment. Regarding the times reported as 0.000: Maple can't measure execution times less than 0.015 seconds; they just get reported as 0. If you need them measured (which doesn't seem to be the case at the moment), then you need to accumulate several runs of the code in one profile.

Next, I decided to measure the effect of precision (as set via Digits) on performance. I set Digits to 30 and redid the above. (Extremely important: This is done starting with restart. Without that, the measurements are difficult or impossible to interpret. I simply don't show the whole thing here because other than Digits, the code is exactly the same.) Here's what I got:

restart:
#omitted procedure code ....
Digits:= 30:
#omitted code that sets up profiling ....
Erweiterung(-1,1,2*x^2+2,2/sqrt(1-x^2),[-.8660254037, 0, .8660254037], 4);

Error, (in Erweiterung) unable to compute coeff

So, you'll need to investigate that error further; and I'll look into it also. Nonetheless, the Profiling information collected up until the point of the error has still been collected, and in this case it's highly relevant: 

Erweiterung
Erweiterung := proc(Unten, Oben, f, g, Liste, n)
local Basenwechsel, p, i, c, d, e, Hn, Koeffizienten, a, s, j, M, V, S, K, nNeu, Em,
Hnm, KnotenHnm, KoeffizientenHnm, h0, b, gxi, Gewichte, Delta, Ergebnis, 
  Endergebnis, Koeffizient, Rest, Basenwechsel1, Basenwechsel2, k;
     |Calls Seconds  Words|
PROC |    1  70.578 932086456|
   1 |    1   0.000     37| Basenwechsel1 := proc( Dividend, m ) :: list;
                                local Basenwechsel;
                                Koeffizient := quo(Dividend,p[m],x);
                                Rest := rem(Dividend,p[m],x);
                                if m = 0 then
                                    Basenwechsel := [Koeffizient*evalf(Int(g*p
                                      [m]^2,x = Unten .. Oben))];
                                else
                                    Basenwechsel := [Koeffizient*evalf(Int(g*p
                                      [m]^2,x = Unten .. Oben)), op(procname(
                                      Rest,m-1))];
                                end if;
                            end proc;
   2 |    1   0.000     31| Basenwechsel2 := proc( Dividend, m ) :: list;
                                local Basenwechsel;
                                Koeffizient := quo(Dividend,p[m],x);
                                Rest := rem(Dividend,p[m],x);
                                if m = 0 then
                                    Basenwechsel := [Koeffizient];
                                else
                                    Basenwechsel := [Koeffizient, op(procname(
                                      Rest,m-1))];
                                end if;
                            end proc;
   3 |    1   0.000    273| p[-1] := 0;
   4 |    1   0.000     10| p[0] := 1;
   5 |    1   0.000      6| for i to 2*numelems(Liste)+2*n do
   6 |   12  70.578 932082225|     p[i] := x^i-add(evalf(Int(x^i*p[j]*g,x = Unten
                                  .. Oben))*p[j]/evalf(Int(p[j]^2*g,x = Unten
                                  .. Oben)),j = 0 .. i-1);
   7 |   12   0.000      0|     userinfo(1,Erweiterung,p[i]);
   8 |   12   0.000    687|     c[i-1] := coeff(p[i],x,i)/coeff(p[i-1],x,i-1);
   9 |   11   0.000    993|     d[i-1] := (coeff(p[i],x,i-1)-c[i-1]*coeff(p[i-\
                                  1],x,i-2))/coeff(p[i-1],x,i-1);
  10 |   11   0.000     33|     if i <> 1 then
  11 |   10   0.000   1886|         e[i-1] := coeff(p[i]-(x*c[i-1]+d[i-1])*p[i
                                      -1],x,i-2)/coeff(p[i-2],x,i-2)
                                else
  12 |    1   0.000    275|         e[i-1] := 0
                                end if
                            end do;

Every line after that has an execution count of 0 due to the error. We see that line 6 is still the time hog. Now look at lines 8 and 9. There is no conditional-statement branching between these lines, yet line 9 has an execution count (11) one less than line ​​​​​​8 (12). How can that be? It's because the error must've occured on line 9.

Like this:

seq([seq](p), p= Iterator:-CartesianProduct([1,2,3] $ 2));

Let's say that you've already created the basic plot and assigned it to the name P:

P:= plot(...);

And let's say that your special point is (x,y), Then do

plots:-display(
    P, 
    plot(
        [[x,y]], 
        style= point, symbol= solidcircle, symbolsize= 20,
        color= red
    )
);

Here's a small modification to acer's final version that uses the very old but little-known type &under:

foo:= proc(
    ode::{`=`, And({set,list}(`=`), Not(0) &under nops)},
    func::function(name), 
    $
)  
    #...
end proc:

type(e, T &under F) is equivalent to type(F(e), T). I'm not posting this because I think this is any better than acer's; rather, I just like to promote the use of this type. "Not zero under nops" is very understandable to me.

By the way, the quotes that you had on `function` do absolutely nothing, neither good nor bad. So, I hope that you weren't expecting them to do something.

The following will get you most of the way to your final form. To finish, you need to apply the well-knowm "half-angle identity" for sqrt(1-cos(x)). Maple has a strong preference for cos instead of sin, and I can't persuade it to change its mind about that.

R:= 1/n*sum(exp(I*(i-1)*k*varepsilon), i= 1..n):
simplify(abs(R)) assuming real, n::posint;
                                           (1/2)
               (-2 cos(k n varepsilon) + 2)     
               ---------------------------------
                                          (1/2) 
               n (2 - 2 cos(k varepsilon))      

The command print never produces programmatically accessible output; the only way to access it is copy and paste.

I can't see any reason why you'd want your output as a set, although it's possible that you have a good reason.

Here are two syntax alternatives for the creation of a column Vector. The first two lines require relatively recent versions of Maple:

f:= x-> local y; fsolve(x*y=1):
R1:= <(for x from 1 by 0.2 to 2 do f(x) od)>;
R2:= <seq(f(x), x= 1..2, 0.2)>;
 

You need to correct the syntax in ode2. After that your dsolve with D(u)(0)=0 works perfectly:

ode2:= -v + mu*diff(r*diff(u(r), r), r)/r = 0:
dsolve({ode2, u(R)=0, D(u)(0)=0});
                                 2      2
                              v r    v R 
                       u(r) = ---- - ----
                              4 mu   4 mu

My corrections are shown in green. And, just to improve readability, I removed three levels of unnecessary parentheses.

It must be noted that while D(u)(0) = 0 implies u(0) is finite, the converse is far from true. D(u)(0) = 0 is a quite strong condition. 

The syntax for a multiple integral requires square brackets for the 2nd argument (the variables and ranges), like either of these:

simplify(Int(f(x,y,z)*Norm(n), [u= 1..2, v= 1..4]));
int(f(x,y,z)*Norm(n), [u= 1..2, v= 1..4], numeric);

However, your surface is undefined along the line y=3/4*x, which unfortunately passes through the rectangle (x= 1..2) x (y= 1..4). By some rearrangements of the variables, you can persuade int to give an answer, but that answer is Float(undefined).
 ​

This seems to me to be a bug in solve. The command eliminate, already shown by acer, is an algebraic generalization of solve, although it has far fewer options.

eqs:= {x*y^2-x, x*sin(Pi*y)}:
map2(eliminate, eqs, {x,y}); #Solve for each variable independently
         {[{x = 0}, {}], [{y = -1}, {}], [{y = 1}, {}]}

#Any return from eliminate whose 2nd element is not the empty set
#is not a complete solution (in the sense of solve). The possible
#presence of these incomplete solutions is what I mean by "algebraic
#generalization". There are no such solutions in this case; I include
#the next line only for the sake of completeness.
(C,InC):= selectremove(s-> op(2,s)={}, %);
      C, InC:= {[{x = 0}, {}], [{y = -1}, {}], [{y = 1}, {}]}, {}

#Remove the superfluous empty sets:
op~(..-2, C);
                  {{x = 0}, {y = -1}, {y = 1}}

#If you want to include the free variable in each solution:
map(solve, %, {x,y});
       {{x = 0, y = y}, {x = x, y = -1}, {x = x, y = 1}}

We must assume that discont returns all of the discontinuities, although it's pretty easy to have one's faith in this shaken. IMO computer algebra is not possible (at this time at least) without tacitly assuming that functions are piecewise continuous. 

The "flavors" used in RandomTools:-Generate can be nested arbitrarily. So, using the Vector flavor,

V:= RandomTools:-Generate('Vector'(variable(length = 1), 4));