The problem is that when you load a Units subpackage (such as with(Units:-Simple)), all of the basic arithmetic operations are replaced with units-aware versions. This is one of the hazards of the with command. The replaced operations are listed as part of the output of the with command. One of those operations is max. I don't know why the Units version doesn't work in your case, but you can call the original max with :-max. So, you can use :-max~(A, B), where A and B are matrices, or you can avoid the need to create a zero Matrix by using :-max~(0, B). The :- prefix can be used on any name to guarantee that it refers to a global version.

Perchance, have you created your own units? I didn't look at the code underlying your embedded components because I don't know how to.

A homogenous linear constant-coefficient recurrence can be done quickly by matrix powering:

((<0,1,0; 0,0,1; 1,-2,1>^(50-3)).<2,-1,0>)[-1];

The last row of the matrix is the coefficients. The other rows are always a staircase pattern of elementary unit vectors as shown. The 3 is the number of coefficients. The vector is the initial conditions. The -1 selects the last entry of the result.

If the exponent is large, say 100,000, then this method is significantly faster than methods that require computing the entire preceding sequence of terms.

A:= proc(n::And(posint, Not(1)))
local a, a1:= 1.75, a2:= 2.25;
    to n-1 do 
        a:= sqrt(a1+4)/6/a2;
        (a1, a2):= (a, a1)
    od;
    a
end proc
:
A(182);

This avoids the need to store the entire sequence of previously computed values.

Use a procedure, like this:

MySum:= (j::integer)->
    add(combinat:-numbcomb(j+i, i+1)*w^(j-1)*(1-w)^(i+1), i= -1..j-2)
:

Like this:

M:= module()
local
    age_type,
    ModuleLoad:= proc()
        TypeTools:-AddType(age_type, integer[0..149])
    end proc
;
export
    CheckAge:= a-> evalb(a::age_type)
;
    ModuleLoad()
end module
:
M:-CheckAge(50);
                              true

type(50, age_type); 
Error, type `age_type` does not exist

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.