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

Please do not substantially change your Questions after they have been answered. You have removed almost all the context that makes Preben's and my Answers understandable. For example, Preben refers to "parts c) and d)". How can a future reader understand that if you've removed those parts from the Question? Thus, you have reduced the value of our work. Is that fair, or ethical? Does that make sense?

Please put back the question as best as you can with your attached worksheet. If you won't, then you're not welcome to ask any more questions here as far as I'm concerned.

@Lottie 

Yes, it works for contour plotting commands. Don't be afraid to try these things yourself even if you're not  sure if they'll work.

Preben has given a detailed analysis of the stability. Here's the quick way, which can be computed automatically with Maple: For nonlinear stability at the origin it is sufficient that the function that we just plotted (dV/dt) have a local maximum at the origin. Obviously this can be seen to be true from the plot. We can verify that computationally. If you recall the multi-variate second-derivative test from Calculus III, for a local maximum, it is sufficient that its gradient be 0 and that its matrix of second partial derivatives (wrt x and y) (called the hessian) be negative definite (i.e., have all negative eigenvalues) at the origin. Thus

v:= [x,y]:
DV:= subs(v(t)=~ v, dV);
VectorCalculus:-Gradient(DV, v);
is(norm(VectorCalculus:-evalVF(%, <0,0>)) = 0);
VectorCalculus:-Hessian(DV, v= [0,0]);
LinearAlgebra:-IsDefinite(%, query= negative_definite);

I just found an amazingly well written set of "notes" on Lyapunov stability analysis. This is better than most textbooks. The author is Gunnar Soderbacka; the link is http://www.users.abo.fi/gsoderba/PhaseP/ljap13.pdf

@acer The Question originally had a link to a journal article---a link which I followed. The article was from a Springer journal dedicated to BVPs. It was dated 2012. The lead author's surname was, I think, Aylin. Amazingly enough, that crap code came verbatim from that article. These authors and their editors have no respect for the code! Obviously it wasn't vetted because it won't even execute. I found the same errors as you did. And there are numerous other issues that are not outright errors for this case, but which could become errors if one tries to adapt the code to other cases.

@student_md (OP): Please put that link back, unless some legal reason prevents you from doing so. I think that since we're seriously criticizing the article, it's considered "Fair Use" under copyright law. (Disclaimer: I have no status or right to offer actual legal opinions.)

@tomleslie I believe that my Answer does what the OP intended rather than duplicating what the OP's code actually does. He did afterall say "The conditions I tried were not working." To me, that means "The inequalities in my piecewise are incorrect."

Furthermore, I believe that he was trying to use piecewise only to do a job that is better done by min, likely because he was unaware of the existence of min.

If I got that wrong, I'll be happy to correct it.

Please clarify: Do you mean that if the return value otherwise computed by the piecewise function is greater than 10, then you want simply 10 returned?

@Rohith You wrote:

• Your first solution is working perfectcly fine almost all cases except when when I have  complex terms in denomenator like in case 4. I would be glad if there is way to solve these.

??? Well what about my revision two Replies above this, "It's not flat"? Doesn't that provide a "way to solve these" for your example 4??  When I post an update or revision to some code, please discard any previous versions and focus all criticism on the the most-recent version. I can't maintain separate versions.

Three Replies above, you said that you wanted the output as 

a*cos(psi)*sin(a)^2, a

---in other words, the function term first, then the coefficient, and with the coefficient attached to the term. But now you show that you want the coefficients first, then the term, and with the coeffcients removed from the terms.

Here is a revision. It handles all five examples from your Reply immediately above. In particular, it handles the distribution of 1/b from example 5.

getFuncCoeff:= proc(e)
option `Author: Carl Love <carl.j.love@gmail.com> 20-Nov-2018`;
local 
     One, sav, FRZ, FRZn:= 0, Freezer:= table(),
     #The thaw of ordinary freeze/thaw is not subtle enough for this.
     Freeze:= proc(e) 
          FRZn:= FRZn+1;
          Freezer[FRZ[FRZn]]:= e;
          FRZ[FRZn]
     end proc,
     #Thaw and recompose the function:
     Thaw:= proc(f::function) 
          local fn:= op(0,f), fc:= Freezer[fn];
          op(1,fc)(op(2,fc)(op(f)))
     end proc,
     Funcs:= [select(
          x-> membertype({function, function^anything}, {op(x)}),
          indets(
               subsindets(
                    subsindets(
                         frontend(expand, [e]), 
                         anyfunc(function),
                         ff-> Freeze([op(0,ff), op([1,0], ff)])(op(op(ff)))
                    ),
                    function, 
                    f-> One*f
               ), 
               `*`
          )
     )[]],
     Coeffs:= subsindets(Funcs, {function, identical(One)}, 1)
;
     Funcs:= eval(Funcs/~Coeffs, One= 1);
     do 
          sav:= Funcs; 
          Funcs:= subsindets(Funcs, 'typefunc'('specindex'(FRZ)), Thaw)
     until sav = Funcs;
     Coeffs, Funcs
end proc:

But you've still not made clear how to handle

1/b*(v1*sin(x)+v2)*(v3*cos(x)+v4).

How should the 1/b be distributed?

Also, you haven't made clear how to handle functions of multiple arguments where some arguments are themselves functions and some other arguments are not. An example would be

arctan(sin(x), 1+cos(x))

When you post examples, please post them in plaintext form so that I can copy-and-paste them.

Your attached worksheet contains the equations only in output form. Ordinarily I'd be able to convert that to input form, but I can't do that with yours. I think there may be a size limit above which that doesn't work and yours is above it. So, could you please supply your equations in input form so that I can work with them?

@666 basha I have something for this situation, but not a complete solution yet. This will work for any 3D surface plot, not just those in my worksheet. If you include the options shading= ZHUE, lightmodel= NONE, you'll get something like this:

ParamPlot3d(
   GetNu, S= 0..12, Ha= 0..12, [Rd= 4.5, Ec= 0.06, Nb= 0.15, Nt= 0.15, Sc= 0.5],
   labels= [S, Ha, Nu],
   surfaceopts= [shading= ZHUE], lightmodel= NONE
);

You will notice that this has not changed the underlying contourplot. You can take any 3D plot and view it as a 2D plot by including the option orientation= [-90,0]. So, for the above

ParamPlot3d(
   GetNu, S= 0..12, Ha= 0..12, [Rd= 4.5, Ec= 0.06, Nb= 0.15, Nt= 0.15, Sc= 0.5],
   labels= [S, Ha, Nu], 
   surfaceopts= [shading= ZHUE], lightmodel= NONE, orientation= [-90,0]
);

I have modified the code so that all dsolve solutions are stored in a remember table; so, it's possible now to change the options on a plot and have the new plot come up instantly. The new worksheet is nanofluid_BVP_Sheikholeslami.mw  The example above is the first of the 3D plots.

Maple's 3D plotting is more-sophisticated than its 2D plotting. Translating the coloring options of a 3D plot to an actual 2D contour plot is difficult (but I think that it's ultimately possible). The built-in coloring option of a contourplot allows for a range bounded by two colors, but not the whole red-to-violet visual spectrum (aka the "HUE" spectrum in some Maple commands). If the above meets your needs, great. If not, I'll try to modify the 2D plot also.