@Lottie The purpose of procedure Q is to defer the evaluation of the expression Q(A) (which is needed in the animate command) until numeric (realcons) values are supplied for A. Those values will be supplied by animate due to its argument A= 2..3. When a numeric value is given, it'll be used to set the parameter A that occurs in the ODE system. Note that there are actually three different variables named A in use here, used a total of seven times:

1. The A in the ODE system and the parameters argument to dsolve;
2. The three occurences of A in procedure Q;
3. The two occurence of A in the animate command.

It's not necessary that the same symbol 'A' be used for these three. But that's a natural way to code it. They do all represent exactly the same thing, but Maple considers them different.

@adel-00 For a parametric plot, the third member of the list that is the first argument of plot must be of the form name=range(realcons). So, you can't have sqrt(Y)= 0..1; you must change it to something like Y= 0..1.

@Carl Love My answer above applies to the general case of an equation in three variables. Since it can be solved for one of the variables, epsilon, as shown by Kitonum below, this is the better way for this case. It's almost always better to solve for a variable than to use implicitplot3d.

@adel-00 The problem is that the contourplot command is not prepared to handle an implicit expression. But implicitplot3d can plot contours, and if these are viewed correctly, they'll look like a 2D contourplot.. So, repeat the same implitcitplot3d, and add the options grid= [49$3], style= contour, orientation= [-90,0]. The orientation makes it look like a 2D plot, but if you pick it up and rotate it, you'll see that it's not.

@adel-00 Please give the expressions that set the  numeric values of the lambdas. I know  that you had those in your separate Question, but that's deleted now.

@mostafaorooji It can be done with coeffs like this:

[coeffs(rhs(ans), indets(ans, suffixed('_C', integer)), 'C')];

This works regardless of the presence of any specific function such as exp. The returned order will correspond to the integers affixed to _C, (which are not necessarily consecutive). If there's any doubt about the order, it will be explicitly listed in C.

A Maple module is a piece of high-level Maple code, essentially a type of procedure. Modules are not stored as machine-language code, and are not stored in .dll files. Some modules may have a .dll file associated with them, but I didn't want you to think that there was an inherent connection between the two.

@vv Huh. I wonder if that's a change made for Maple 2018.2. I still have 2018.1, and RealRange(-infinity, infinity) returns unevaluated. Anyway, I changed the code to accomodate this situation, so please test it.

@vv Please apply lprint to the results containing real, and tell me what you get.

@Harry Garst It would not have occurred to me to use this method as a starting point for factoring.

@adel-00 It's just a trivial change:

plot([rhs(f)^2, Y, Y= 0..10], labels= [typeset(epsilon^2/`4`/Pi), Y]);

@vv I corrected the code to accomodate your ex1, and the new code is in my previous Reply. Regarding your ex2, the result that I get is RealRange(-infinity, 1). Unless I'm having a huge mental lapse, that's the correct answer. Is that what you get?

@acer I did both things, and here's the final product:

SimplifyIntervals:= module()
description
      "Simplifies arbitrarily nested set expressions of real intervals expressed with RealRange"
      " and/or finite discrete sets. Returns a union of disjoint intervals and discrete points."
;
option `Author: Carl Love <carl.j.love@gmail.com> 26-Nov-2018`;
local
       x::identical(x), #bound variable for `solve`
       
       #1. Declare a type specific to this module.  2. Hook into `simplify`.  
       ModuleLoad:= proc($)
          TypeTools:-AddType( #recursive type: put base cases 1st!
                 'RealIntervals', 
                 Or(  
                    specfunc(RealRange), set(realcons), identical(real), #base cases 
                   'specop'('RealIntervals', {`intersect`, `union`, `minus`}) #recursive cases
                )
          );
          :-`simplify/intervals`:= ()-> ModuleApply(args);
          forget(simplify)
     end proc,
     
     #Convert set expression into boolean expression with bound variable x.
     Deconstruct:= (IN::RealIntervals)->
          if IN::specfunc(RealRange) or IN=real then convert(x::IN, relation)
          elif IN::set(realcons) then Or((x =~ IN)[])
          elif IN::`minus` then And(thisproc(op(1,IN)), Not(thisproc(op(2,IN))))
          else `if`(IN::`union`, Or, And)(map(thisproc, [op(IN)])[]) 
          fi,

     ModuleApply:= proc(IN, {no_union::truefalse:= false})
     #Using 'no_union' option guarantees that output will be sorted the standard way, even if it contains 
     #discrete points.
     local P, J;
          if not IN::'RealIntervals' then return IN fi;
          #Separate 'solve's output into discrete points (P) and 'RealRange's:
          (P,J):= selectremove(type, [solve(Deconstruct(IN), x)], realcons);
          if x in J then return 'RealRange'(-infinity, infinity) fi;
          J:= sort([`{}`~(P)[], J[]], 'key'= (p-> `if`(op(1,p)::realcons, op(1,p), op([1,1],p))));
          `if`(no_union, J, `union`(J[]))      
     end proc
;
       ModuleLoad()
end module: 

Example (showing three formats for the output):

IN:= (((RealRange(-infinity, 1) intersect RealRange(Open(-1), infinity))
   minus RealRange(-1/3, 1/3)) minus {0}) union {-2,2};

simplify(IN, 'intervals');
SimplifyIntervals(IN, no_union);
lprint(%);
 

@acer Thanks. I did exactly what you did above, except that I forgot to forget. 

Thank you for at least restoring some parts of the Question.