@vv  Thanks for this method. The key to success here, of course, is to rewrite the system as a polynomial system (without  sqrt ). Unfortunately, we have to do it manually. After that, the  solve  command also solves the problem.

sys := r-c>0, 0<=r^2 - 2*r*c -r +2*c, r^2 - 2*r*c -r +2*c  < (r-c)^2, 0<=c, c<=1 , 0<r, r<1, r<2*c:
solve([sys], r, parametric=true);

Wolfram Mathematica solves the problem immediately in its original setting:

@manju  1. Maple can choose a random object only from a certain finite set. Therefore, a certain range must be specified.
2. Replace colons with semicolons, for example  a_1:=rand(1..5)();  and so on.

It is enough to use backticks only in denominators, for example   4=456/`1920`  and so on.

@Maple_lover1  I corrected my answer.

@BrettKnoss  See my solution above.

@BrettKnoss  I didn't understand the meaning of your question. If a series is given, then the general term of the series  u(n)  must be known for any  n .

@BrettKnoss  c(n)  is just the coefficient in front of  n_th  term of the series:

restart;
S := Sum((x-2)^n/(n^2-1)^2, n = 2 .. infinity);
op(1,S);
coeff(%, (x-2)^n);

@matviiv10 

restart;
h:=0.01: r:=0.02: N:=10: f:=y->y^5: K:=(x,y)->sin(x)*y^5:
for n from 1 to N do 
x[n]:=sin(h*n);
od:
for m from 1 to N do
y[m]:=cos(r*m);
od:
n:='n': m:='m':
Sys:={seq(sum(K(x[n],y[m])*U[n],n=1..N)=f(y[m]), m=1..N)};
solve(Sys);
assign(%);
plot([seq([x[n],U[n]], n=1..N)]);

@ahmeng  Above, I showed only one solution (in my opinion, the most direct and reliable). Here are 2 more ways:

Download more_ways.mw

@ibrahimcbx 

1. You can do this with either separate commands or with a single command.

2. I didn't understand the meaning of your second question.

@Carl Love  It seems that your solution will be correct only in the segment  p in [0, 3-2*sqrt(2)]

@salma951  I didn't notice that 3 appears twice and 6 is missing. So write  randperm([1,2,3,3,4,5]) 

@jalal  Replace  range = -3 .. 10  by  (for example)  range = 1 .. 10  .

@acer  I think in all intermediate calculations it is better to use exact mathematics. Of course, we can present the final results in any desired form.

For example:

[solve({eqy3, x>=-5, x<=5})];

nm 4588  We can easily generalize your method to other exponential functions (in addition to exp). We can also use the  applyop  command (many shorter if the position of the corresponding terms is known):

restart;
expr := -2*exp(x)*x - exp(x) + 2*exp(x)^2*x + exp(x)^2 +(3^x)^2+ sin(x)^3 + cos(x)^2;
subsindets(expr,{exp(anything)^anything,(anything^anything)^anything},combine); # The first way
applyop(combine,{$1..5},expr); # The second way