@binbagsss You wrote:

•  I will need to define an alias

One never needs to define an alias. You can use direct assignment instead, as long as you don't use the same names on both left and right sides of the assignment operator. The only potential benefit of an alias (in this situation) is the way that it can simplify the output. However, I think that that simplified output is causing you some confusion. So, I think that for the time being you shouldn't use alias.

An example of potentially confusing aliased output is D(V). You've aliased V to V(h(x)). So does D(V) then mean D(V(h(x)))? An experienced user probably realizes that this latter expression is pretty much nonsense, and thus knows that that is not what was meant. Thus they realize that the V in D(V) is just plain, original V. 

So what does D(V) mean? It considers V to be a univariate functional operator and is expressing the operator form of V's derivative with respect to its unnamed variable. Some more-explicit examples should help: 

D(sin) returns cos.
D(sin(x)) is nonsense, because sin(x) is an expression, not an operator.
D(sin)(x) returns cos(x).
D(sin)(0) returns 1.
diff(sin, x) is nonsense, because sin is an operator, not an expression depending on x.
diff(sin(x), x) returns cos(x).
diff(sin(1), x) returns 0 because sin(1) is constant with respect to x.
diff(sin(1), 1) is an error, because 1 is not a variable.
D(sin)(1) returns cos(1).

The error message that you show is particularly common with boundary-layer BVPs. Unlike most error messages, it is not a result of a programming error, but represents something inherent to the structure of the BVP. Often it means that there are multiple solutions. You have parameters br, k1, and lambda. (You also have parameter blt as a "fake infinity", but, for the moment, I don't want to consider that to be a changeable parameter.) The way around the error is to use approxsoln, as mentioned by @dharr . The first step towards doing that is to experiment with varying the parameters (to "reasonable"^[*1] values) until you find values for them for which you don't get an error. Then we can gradually and contuously shift them to the parameter values that you want with each step using the previous step's result as its approxsoln. If you can find the values for the first step, I'll help with the rest.

[*1] Setting all parameters to 0 would likely converge without error, but that's probably not reasonable.

The key to an efficient solution procedure is to solve the equation symbolically (with solve) just once, and then construct all subsequent numeric solutions from that. Like this:

eq:= %abs[i](%a*x + %b) + %abs[j](%c*x + %d) - %t*x^2 + %m*x - %n = 0;

X:= subsindets(
    {seq}(
        seq(solve(eval(eq, %abs= [x-> x, x-> -x]), x), i= 1..2), 
        j= 1..2
    ),
    anything^(1/2),
    s-> %sqrt(op(1,s))
);

soln2:= subs(
    _X= X, [%a, %b, %c, %d, %t, %m, %n]=~ [a, b, c, d, t, m, n],
    (
        a::integer, b::integer, c::integer, d::integer,
        t::integer, m::integer, n::integer
    )-> 
        (X-> if nops(X)=6 then [[_rest], X] else fi)(
            select(type, value(_X), integer), args
        )
):

Sols:= CodeTools:-Usage((soln2@op)~(convert(L, listlist)));
memory used=29.38MiB, alloc change=-4.00MiB, 
cpu time=360.00ms, real time=334.00ms, gc time=78.12ms

     [[[5, 4, 3, 7, 1, 1, 1], {-4, -3, -2, -1, 1, 10}], 
         ...40 other solutions elided... 
     ]

I'm curious about something: Your equation can obviously have up to 8 solutions. Why are you only interested in 6?

Here is an Answer to your explicit Question about map: The command map replaces only a single function argument (although not necessarily the first argument) with values from a container (vector, list, set, etc.). To do the equivalent thing with multiple arguments simultaneously, use the elementwise operator ~. In your particular case, you could do

(tgf@op)~(convert(L, listlist))

where L is a Matrix or 2D Array. The following also do it:

(tgf@seq)~(convert(L, listlist))
map(tgf@op, convert(L, listlist))

and many other variations are also possible.

However, if your primary concern is reducing time, a change along these lines is insignificant. The time for operating the looping mechanism (whether it be for, seq, map, ~, etc.) too small to even be measured; the vast, vast majority of the time is used by solve. But you don't need to use solve. I gather from your worksheet that you're interested in finding integer 7-tuples (a, b, c, d, t, m, n) for which the equation abs(a*x + b) + abs(c*x + d) - t*x^2 + m*x - n = 0 has 6 integer solutions for x and then listing those solutions. Is that correct? If so, it can be done much more easily than by using solve, fsolve, or any other real- or complex-valued root-finding technique. I'll do it in a separate Answer.

First, a note on terminology: An integral equation is an equation that contains an unknown function inside an integral. What you have is just an integral, which is much simpler than an integral equation.

You have used Sum instead of sum, which expresses the sum in inert form. That is, the sum is shown in Sigma notation (with large uppercase Greek sigma), even though in this case it's trivial to expand the sum to its terms. There's no problem with using an inert form; indeed, it makes the displayed expression much more readable. The presence of inert forms is indicated by a gray (rather than blue) character (or characters) in the prettyprinted display. In this case, the Sigma is gray. Likewise, Int can be used as an inert form of int.^[*1] That wouldn't make much difference in this case. 

Kitonum's Answer replaces Sum with add. This can be done (and often should be done) in any case where the sum has a specific finite number of terms. However, you lose the benefit of the prettyprinted display. You also lose the option of the sum being simplified symbolically via a summation rule (that doesn't apply in this specific case).

The command value can be used to convert inert forms to their non-inert counterparts. So, in this case, you'd use value(eq1).

[*1] Don't assume that every command has an inert form obtained by capitalizing it. The commands for which this is true are diff, eval, int, limit, product, and sum. An inert form of any command foo (including those already listed) is %foo.

Try this link: Maple Learn

The GUI (user interface) is completely different than any of those offered by regular Maple, but I think that the "kernel" (mathematical engine) is the same.

Like this:

restart:
f1:= (x,y,z)-> x*y*z:
f2:= (x,y)-> x+y:
(x,y,z):= (2,3,4):
Grid:-Run(0, f1, [x,y,z], set= {f1}, assignto= 'r1'):
Grid:-Run(1, f2, [x,y], set= {f2}, assignto= 'r2', wait):

r1, r2;
                             24, 5

As you probably realize, the wait is not necessary if you have other work that could be performed asynchronously before r2 becomes available.

Use m &^ e mod n.

The reason that the problem occurs is the order of operations: In m^e mod n, the m^e is computed (or at least an attempt to compute it is made) before that number is passed to mod. But m &^ e is returned as an unevaluated function call `&^`(m, e), and that is exactly what gets passed to mod.

Here's my rewrite of GraphTheory:-UnderlyingGraph:

Underlying:= proc(G::GRAPHLN) 
uses GT= GraphTheory;
local
    V:= GT:-Vertices(G),
    Vi:= table(V =~ [$1..nops(V)]), 
    W:= subs(
        _W= (W-> W+W^+)(GT:-WeightMatrix(G)), 
        (u,v)-> `if`(nargs=1, _W[Vi[u]$2]/2, _W[Vi[u],Vi[v]])
    ),
    e
;
    GT:-Graph(V, {seq}([e, W(e[])], e= (`{}`@op)~(GT:-Edges(G))))
end proc
:

In the interest of time, I've omitted some of the bells & whistles of the original procedure, focusing instead on precisely the features that you implied that you needed.

I am investigating the bug in the original procedure.

Your worksheet shows a 1st-order ODE. Thus dsolve will only use one boundary condition. You said that A is a constant of integration. Well, where is the original 2nd-order ODE that produced that constant? That's what you should pass to dsolve, along with both boundary conditions.

You wrote:

• [I]t seems like the numerical solver cannot solve this type of nonlinear ODE.

It's not that it cannot solve it; it's simply that it chooses not to solve it unless you provide further information that clarifies an ambiguity. Whichever branch of the ambiguity you choose, the numerical solver can very easily handle it. Your ODE is

y' + (y')^2 = 0.

Trivial algebra (which can be done with the solve command if you wish) shows that that ODE is equivalent to

y' = 0  or  y' = -1.

You need to pick one of those choices before passing it to dsolve(..., numeric).

You'll get the same error in any case where the ODE cannot be uniquely solved (algebraically, not as a differential equation) for its highest-order derivative. As long as that can be done, the numeric solvers make no distinction between linear and nonlinear ODEs. However, since your stated interest is nonlinear, it should be pointed out that in this particular, trivial case, once you clarify the ambiguity, the ODE becomes linear.

In a Reply to @nm above, I mentioned the need for the student to begin with "geometric intuition" for this general class of problem. Both the Answer by Kitonum and the one by acer show great use of this intuition. However, in cases where Green's Theorem can be applied, very little intuition is needed. The Theorem can be applied when the curves can be parameterized; then the area integral can be replaced by a line integral. In other words, only the boundary matters (which is rather amazing), no matter how reticulated it is. (It's as if I could measure the area of the pond in front of my house just by walking around it's boundary, never looking across the water).

The limits of integration are determined by the intersection points alone; there's no need to find roots of the cubic.

In the following transcription, I've deleted  much of the output consisting of very long algebraic numbers (some of very high degree). You can see it in the worksheet. The given solution is exact.

restart:
Digits:= 15:
eq1:= x^2/5 - y^2/2 = 1:  eq2:= y = x^3 - 7/2*x^2 + 2*x:

#Parameterizations (must use a different parameter for each):
C1:= table([x= sec(s)*sqrt(5), y= tan(s)*sqrt(2)]):
C2:= table([x= t, y= subs(x= t, rhs(eq2))]):

#Find intersection points:
solve({C1[x]=C2[x], C1[y]=C2[y]}, {s,t});

#Filter out any that are not real:
R:= select(type, {allvalues(%)}, set(name=realcons));
#Count the intersections:
nops(%);
                               2

#Great! Having only 2 intersections makes this especially easy!

#We don't need float values for the rest of this problem, but just out
#of curiousity...
Rf:= evalf(R)
    Rf := {{s = -0.540880221771561, t = 2.60840228710112}, 
      {s = 0.716707677716986, t = 2.96571533700056}}

#We can now plot directly rather than using implicitplot:
plot(
    [
        [C1[x], C1[y], s= eval(s,R[1])..eval(s,R[2])],
        [C2[x], C2[y], t= eval(t,R[1])..eval(t,R[2])]
    ],
    thickness= 3
);

#The plot shows a simple connected region. Good!

#The following integral will either give the area or its negative,
#I don't care which.
A:= 
    int(C1[x]*diff(C1[y],s), s= eval(s, R[1]) .. eval(s, R[2])) +
    int(C2[x]*diff(C2[y],t), t= eval(t, R[2]) .. eval(t, R[1]))
;
evalf(A);
                        0.75833670589040

#I think that A is in an extension field of the algebraic numbers generated by 
#including the ln of a few of them, but I'm not sure.

op~(0, indets(A, function));
                          {RootOf, ln}
#A is extremely big:
length(A);
                             52898

GreensThm.mw (Download worksheet)

P.S.: My plot above taken by itself doesn't exclude the possibility that one of the curves doesn't re-enter the region for some unshown values of s or t, which would substantially complicate the problem. A broader plot is useful as acer suggests.

Please read carefully the section "Using Assumptions" in the help page ?solve,details.

Conclusion: Don't use assuming real with solve.

After what acer taught you yesterday, I'd've thought that you'd stop using all-variable assumptions such as assuming real for any command. 

This code uses the extremely fast package LinearAlgebra:-Modular. It is much faster than the compiled code, finishing the 2^25 case in under 2 minutes.

[Edit: The complete corrected code is in a Reply below. It does the 2^25 case in 0.2 to 0.3 seconds.]

BooleMobiusTransform:= proc(V::Vector(datatype= integer[8]))
uses LAM= LinearAlgebra:-Modular; 
local 
    nv:= numelems(V), n:= ilog2(nv), im:= Scale2(1, n-1), istep:= im, 
    jm:= 1, h:= Scale2(im, 1), istart, 
    R:= LAM:-Create(2, h, 0, integer[8]);
    LAM:-Copy(2, V, 1..-1, R, 1..nv);
    to n do 
        istart:= 1; 
        to jm do
            LAM:-AddMultiple(
                2, R, istart..im-1+istart, 
                R, istart+istep..im-1+istart+istep
            );   
            istart+= h
        od;
        im:= Scale2(im,-1); istep:= Scale2(istep,-1); 
        jm:= Scale2(jm,1); h:= Scale2(h,-1)
    od; 
    R
end proc
:
for n from 11 to 19 do
    V:= Vector(2^n, rand(0..1), datatype= integer[8]):
    lprint('n'=n);
    V1:= CodeTools:-Usage(BooleMobiusTransform(V))
od:
n = 11
memory used=0.58MiB, alloc change=0 bytes, 
cpu time=15.00ms, real time=9.00ms, gc time=0ns
n = 12
memory used=1.16MiB, alloc change=0 bytes, 
cpu time=16.00ms, real time=14.00ms, gc time=0ns
n = 13
memory used=2.35MiB, alloc change=0 bytes, 
cpu time=31.00ms, real time=21.00ms, gc time=0ns
n = 14
memory used=4.76MiB, alloc change=0 bytes, 
cpu time=47.00ms, real time=44.00ms, gc time=0ns
n = 15
memory used=9.63MiB, alloc change=0 bytes, 
cpu time=93.00ms, real time=85.00ms, gc time=0ns
n = 16
memory used=19.51MiB, alloc change=-1.00GiB,
cpu time=390.00ms, real time=358.00ms, gc time=93.75ms
n = 17
memory used=39.51MiB, alloc change=1.00MiB,
cpu time=437.00ms, real time=393.00ms, gc time=109.38ms
n = 18
memory used=80.01MiB, alloc change=2.00MiB, 
cpu time=797.00ms, real time=742.00ms, gc time=125.00ms
n = 19
memory used=162.01MiB, alloc change=3.00MiB, 
cpu time=1.62s, real time=1.50s, gc time=265.62ms

Often (but not always), Maple can symbolically integrate the Heaviside forms of piecewise expresssions that it can't do directly. Two branches of your proposed result needed small corrections (assuming I and Maple did everything correctly here).

restart:
interface(showassumed= 0):
assume(n::posint, j::nonnegint, j <= n, T > 0, alpha > 0);
h:= T/n:
psi:= t->
    local i:= op(procname);
    convert(
        if i=0 then piecewise(t<0, 0, t<h, 1-t/h)
        elif i=n then piecewise(t < T-h, 0, 1-(T-t)/h)
        else 
            piecewise(
                t < (i-1)*h, 0, 
                t < i*h,     t/h-(i-1), 
                t < (i+1)*h, (i+1)-t/h
            )
        fi,
        Heaviside
    )
:    
g:= proc()
local i,j,t;
    (i,j):= op(procname);
    simplify(
        int((j*h-t)^(alpha-1)*psi[i](t), t= 0..j*h)/GAMMA(alpha)
    )
end proc
:
g[0,j]() assuming additionally, j > 0;
             1         //       alpha            (alpha + 1)
    - ---------------- \\(j - 1)      (j - 1) + j           
      GAMMA(alpha + 2)                                      

                       alpha\  (-alpha)  alpha\
       + (-alpha - 1) j     / n         T     /


#Compare with your proposed result:
#Note that I added "-" in front of (j-1)!
simplify(
    % - 
    h^alpha/GAMMA(alpha+2)*
        (-(j-1)^(alpha+1) + j^alpha*(1-j+alpha))
);
                               0

g[j,j]() assuming additionally, j > 0;
                         (-alpha)  alpha
                        n         T     
                        ----------------
                        GAMMA(alpha + 2)

#Clearly, that's your proposed result.

g[i,j]() assuming i::posint, additionally, i > j, 
    additionally, i <= n;
                               0

g[i,j]() assuming i::posint, additionally, i < j;
       1         / alpha /           alpha            
---------------- \T      \(j - i - 1)      (i - j + 1)
GAMMA(alpha + 2)                                      

                              alpha            alpha            \ 
   + (-i + j + 1) (-i + j + 1)      + 2 (j - i)      (alpha + 1)/ 

   (-alpha)\
  n        /


#Note how I corrected your formula below:

simplify(
    % - 
    h^alpha/GAMMA(alpha+2)*
        ((j-i+1)^(alpha+1) + 2*(j-i)^alpha*(alpha+1) - (j-i-1)^(alpha+1))
);
                               0