Bisection is the crudest root-finding method that has any practical value, that value being that it's guaranteed to work under some fairly broad conditions (Intermediate Value Theorem) and it's easy to understand. I don't recommend it for this problem, but since you want it, I'll use it anyway. No while loop is needed: Since the interval is reduced by a constant factor on each iteration, the exact number of iterations required can be precomputed. That number is ceil(log[2]((l2-l1)/width)). 

Bisection:= proc(f, int::range(realcons), err::positive:= 10^(1-Digits))
local lo, hi, mid, f_lo, f_mid;
    (lo, hi):= op(evalf(int));
    if (f_lo:= evalf(f(lo)))*evalf(f(hi)) > 0 then
        error 
            "function evaluations at interval ends must be on "
            "opposite sides of 0"
    fi; 
    to 1+ilog2((hi-lo)/err) do
        if f_lo*(f_mid:= evalf(f((mid:= (lo+hi)/2)))) > 0 then lo:= mid
        elif f_mid = 0 then return mid
        else hi:= mid
        fi
    od;
    lo..hi
end proc
:
sol:= dsolve(
    {-p*diff(y(x),x$2) + q*diff(y(x),x) = L*y(x), y(0)=0, D(y)(0)=1},
    numeric, parameters= [p, q, L]
):
solut:= proc(L::numeric)
    sol(parameters= [1, 1, L]);
    eval(y(x), sol(2))
end proc
:
Bisection(solut, 1..5, 1e-6);
                   2.717401504 .. 2.717402458

You need a semicolon (;) after your last end if.

Your ormap command is functional, but much too elaborate. It can be replaced by simply has(sols, _Z). 

The results that you want are returned by evalc(expr), evalc(-expr), evala(expr), evala(-expr), radnormal(expr), and radnormal(-expr).

I guess that there's no straightforward algorithm for general simplification. It might proceed as you would "by hand": by recognizing patterns or tricks in the expression.

There is one method that I know whereby an ODE BVP can be solved using an ODE IVP solution method: the shooting method. Occasionally---often for boundary-layer BVPs such as yours---it gives better resuts than Maple's stock mesh-based BVP solver. Here are the details. There is no need to use 4th-order Runge-Kutta for this; any IVP solver (that accepts parameters) will work.

MaplePrimes refuses to show this worksheet, so I'll transcribe the essential part of the code. The code, comments, output, and plots are in the attached worksheet.

restart;
N1:=1:  N2:= 1:  N3:= 0.1:  R:= -1:
EQ:= {
    (1+N1)*diff(f(x),x$4) - N1*diff(g(x),x$2) -
    R*(-diff(f(x),x)*diff(f(x),x$2) + f(x)*diff(f(x),x$3))
    = 0,
    N2*diff(g(x),x$2) + N1*(diff(f(x),x$2) - 2*g(x)) - 
    N3*R*(f(x)*diff(g(x),x) - diff(f(x),x)*g(x))
    = 0
}:
#Boundary points (lower and upper):
Lb:= 0:  Ub:= 1:

#Separate the boundaries:
BCL:= {(D@@2)(f)(Lb)=0, f(Lb)=0, g(Lb)=0}: 
BCU:= [D(f)(Ub)=0, f(Ub)=1, g(Ub)=0]: #Note: list, not set

#Replace one set of BCs with same number of parametric ICs at the other boundary:
ICLpara:= [D(f)(Lb)=pD1f, (D@@3)(f)(Lb)=pD3f, D(g)(Lb)=pD1g]: #list, not set

dsolIVP:= dsolve(
    EQ union BCL union {ICLpara[]}, numeric, parameters= rhs~(ICLpara),
    method= classical[rk4] #could use any IVP method
);

(BCU_names, BCU_vals):= (lhs~,rhs~)(BCUdiff);

Match:= proc(pD1f, pD3f, pD1g)
option remember;
    dsolIVP('parameters'= [pD1f, pD3f, pD1g]);
    (eval(BCU_names, dsolIVP(Ub)) -~ BCU_vals) #residuals
end proc
:
Match1:= (pD1f, pD3f, pD1g)-> Match(args)[1]:
Match2:= (pD1f, pD3f, pD1g)-> Match(args)[2]:
Match3:= (pD1f, pD3f, pD1g)-> Match(args)[3]:

#This is the iterative part of the shooting method:
fsolve(
    [Match||(1..3)],
    [-9..9, -9..9, -9..9] #search ranges for parameters
);
             [1.5132303350250941, -3.1484135938501150, -0.41182804973022998]

plots:-odeplot(dsolIVP, [[x,f(x)], [x,g(x)]], x= 0..1, axes= boxed);

#Compare with stock BVP solution:
dsolBVP:= dsolve(EQ union BCL union {BCU[]}, numeric);
plots:-odeplot(dsolBVP, [[x,f(x)], [x,g(x)]], x= 0..1, axes= boxed);

Download Shooting.mw

The syntax of the exp function is exp(x) (where x could be any algebraic expression). You have exp^(x).

From the 20 variables, you will need to pick 13 to solve for.

Download PointPlaneDistance.mw

sol:= dsolve(
    {-p*diff(y(x),x$2) + q*diff(y(x),x) = L*y(x), y(0)=0, D(y)(0)=1},
    numeric, parameters= [p, q, L]
):
solut:= proc(L::numeric)
    sol(parameters= [1, 1, L]);
    eval(y(x), sol(2))
end proc
:
plot(solut, 0..100);

fsolve(solut, 1..5);
             2.717402017

@Rouben Rostamian  Here's a simplification of Rouben's procedure:

ZigZag:= (i,j)->
local d:= i+j-1, n:= op(procname);
    `if`(d>n, thisproc(n+1-~(j,i)) + (d-n)*(3*n-d), (d^2+1 + (-1)^d*(i-j))/2)
:
#Example of use:
Matrix(10$2, ZigZag[10]);

Your error message shows that you've used square brackets in the equations. Square brackets cannot be used as substitutes for parentheses.

The command is map(limit, M, [x1= 0, x2= 0]), not the map2 command that you gave. You must respect the order of the arguments when using map and its relatives.

I disagree with Rouben's method. I think that the coeffs command is needed for this:

V:= [u[2], u[2,2], u[2,2,2]]:
C:= coeffs(collect(f, V, distributed), V, 'terms'):
Coeffs:= table(sparse, [terms]=~ [C]);

The following should be faster if you're just as likely to be checking residues as nonresidues:

if NumberTheory:-QuadraticResidue(x,n) then 
    a:= NumberTheory:-ModularSquareRoot(x,n);
    b:= b+1 #or b++ in Maple 2019+
else
    c:= c+1 #or c++ in Maple 2019+
fi;

And if you want the last returned value rather than the last computed value and also use Joe's "last" trick, then do this:

Foo:= ()-> (foo("last"):= foo(args)): #Wrapper procedure
foo:= proc(x) option remember; x^2 end proc:
foo("last"):= "foo not yet used":
Foo(2);
                               4
foo("last");
                               4

Foo(3);
                               9

foo("last");
                               9

Foo(2);
                               4

foo("last");
                               4

Of course you can't apply 1/z or Joukowsky to regions containing z=0.

Here's an example of Joukowsky applied to a triangle:

restart:
PolygonAsLines:= (P::list([numeric,numeric]))->
local t, k;
    plot(
        [seq([((P[k+1]-~P[k])*~t +~ P[k])[], t= 0..1], k= 1..nops(P)-1)],
        _rest
    )
:
J:= z-> z + 1/z; #Joukowsky
Jw:= evalc([Re,Im](J(x + y*I)));
commonopts:= view= [1..2.5, 0..1], size= [500$2], scaling= constrained:
R:= plots:-display(
    PolygonAsLines([[1,0],[2,0],[1,1],[1,0]], color= blue, thickness= 3),
    seq(plot(i, x= 1..2-i, color= black, thickness= 0), i= 0..1, 0.1),
    plot([seq([i, t, t= 0..2-i], i= 1..2, 0.1)], color= black, thickness= 0),
    commonopts
);
F:= plottools:-transform(unapply(Jw, (x,y))):
plots:-display(F(R), commonopts);