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);

By far the easiest way to use parallel programming is through the Seq and Map commands of packages Grid or Threads. You simply take a seq or map command from your sequential code and replace it with the package command:

foo:= i-> i^2:
Threads:-Seq(foo(i), i= 1..9);
Grid:-Seq(foo(i), i= 1..9);
X:= [$1..9]:
Threads:-Map(foo, X);
Grid:-Map(foo, X);

So the complexity of Threads:-Task is not needed here, whereas it may be needed to parallelize code that can't be expressed with seq or map such as recursive code.

Code that writes to global variables (without locking them) is not "thread safe". The vast majority of Maple's high-level symbolic commands do this, thus separate invocations are not truly "independent" (as you said), thus Threads can't safely be used. If you try to use Threads with code that isn't threadsafe, the results are very unpredictable. You may get incorrect results, an error, or a crash.

On the other hand, usage of global variables doesn't affect Grid, but you pay for that with a much higher startup time and memory usage.

@nm z^(-1) obviously has a z in it!  

It helps to think of breaking expressions down into trees:

z^(-1):   Root node:  ^
2 branch nodes:  z  and  -1

There is a z node.

1/(x*y):   Root node:       *
2 branch nodes: x^(-1) and y^(-1)

There's no x*y node.

You could substitute x = t/y.

The command ToInert(1/(x*y)) ​​​​​​or ToInert(1/z) will show you the trees (in 1D function form).