@torabi I made the following naive attempt not thinking of the noncommutativity of differential operators and multiplication operators. See note at bottom. I don't believe it was a problem in your previous example.

I removed the global declaration in px entirely. It is not necessary anyway.

indets(B,name); #We find (besides x,y,z, and pi) X, XY, XZ, Y, Z, ZY.

px := proc (u1) local u;
  u := expand(u1);
  if type(u, `+`) then return map(procname, u) end if;
  if not type(u, `*`) then return u end if;
  if member(X, {op(u)}) then diff(u/X, x) else u end if
end proc;
py := subs(X = Y, x = y, eval(px));
pz := subs(X = Z, x = z, eval(px));
pxx:=subs(X=XX,x=(x,x),eval(px));
pyy:=subs(X=YY,x=(y,y),eval(px));
pzz:=subs(X=ZZ,x=(z,z),eval(px));
pxy:=subs(X=XY,x=(x,y),eval(px));
pxz:=subs(X=XZ,x=(x,z),eval(px));
pzy:=subs(X=ZY,x=(y,z),eval(px));
P:=`@`(px,py,pz,pxx,pyy,pzz,pxy,pxz,pzy);
R:=map(P,B);
###############
### Time for a very critical warning: Differential operators don't commute with multiplication operators.
In your recent worksheet Q has e.g. this element:
Q[2,1];
## In your input you write it as
2/(Pi*a*(a*y+b))*ZY;
## Ordinary multiplication in Maple is commutative, so that could be a problem depending on what your actual meaning was.
This operator will be handled by pzy as
ZY*2/(Pi*a*(a*y+b)); #In that order and was that intended?
## Now compare
diff(2/(Pi*a*(a*y+b))*f(x,y,z),z,y);
## with
2/(Pi*a*(a*y+b))*diff(f(x,y,z),z,y);
## My conclusion is that you must rewrite Q so that the noncommutativity is respected. Maybe use &* or . (matrix multiplication).
The whole must be thought over and another approach used. The problem is not that Q and N are matrices.
You must first be able to handle an operator like 2/(Pi*a*(a*y+b))&*ZY, where ZY stands for diff(..., y,z). I only used &* here to keep the order if entered in Maple.

@mathematicianS In your original formulation of the boundary conditions you have f'''(0), yet the ode for f is only of order 3. You must reformulate that so that only orders lower than the highest derivative order in the ode appears. 
In your first order systems version the same problem appears since D(v)(0) appears, i.e. v'(0).

Closely related questions have been asked and answered often in the past. For anybody to help you give us the odes and the boundary values as text so we don't have to type.
Since solution will be numerical we shall need the values of lambda, Pr, s, sigma, a, and b.

@vv In a programmatic context it appears that the repeated creation of the same procedure keeps the address.
In these examples I removed option remember since it seems to be irrelevant for my point:

restart;
F:=proc() 0 end proc;
addressof(eval(F));
F:=proc() 0 end proc;
addressof(eval(F)); #Different from above
to 10 do F:=proc() 0 end proc; print(addressof(eval(F))) end do: #All addresses are the same
to 10 do F:=proc() 0 end proc; print(addressof(eval(F))) end do: #All the same, but different from just above
F:=proc() 0 end proc; addressof(eval(F)); #Try executing this very line twice. Result: Different addresses.
p:=proc() global f; f:=proc() 0 end proc; NULL end proc; #global might as well have been local, no difference
p();
addressof(eval(f));
p();
addressof(eval(f)); #Same as just above
############
These examples remind me of an unrelated discussion of a weird bug that only occurred when the user entered the code. I have forgotten any details.

@Carl Love Rather strange.
I made this observation. Compare the following 2 procedures.
q1 has your f, q2 has a modified version of your f. f in q2 just has an evaluation of a variable k local to q2.
q1 apparently needs forget, q2 doesn't.
Note added: If the NULL output is changed to f, then it appears that the two successive outputs from q1() evaluate to the same procedure, whereas this is not the case for q2.
Conclusion: The behavior of q1 is not a bug. The procedure with the local name f (name new at its invocation) is the same in each call. The memory table belongs to the procedure, not to its (many) names. Thus the memory table remains.
But the procedure f produced by q2 at each call is so to speak "parametrized" by the local variable k. Thus these procedures f are all different.
####
restart;
q1:=proc() local f, k; #k not used anywhere though
   f:=proc() option remember; 0 end proc;
   print("Entering:", op(4,eval(f)));
   f(89);
   print("Leaving:",op(4,eval(f)));
   NULL
end proc;
q1();
q1();
q2:=proc() local f,k;
   f:=proc() option remember; k; 0 end proc;
   print("Entering:", op(4,eval(f)));
   f(89);
   print("Leaving:",op(4,eval(f)));
   NULL
end proc;
q2();
q2();

@one man Your comment "Unclear" is ironically unclear to me.
If you give us a concrete problem, we can have a look at it.

@Thomas Richard If I do (in Maple 2016.1)

restart;
verify(cos(u)+sin(u), sqrt(2)*cos(u-Pi/4), equal);
## I get FAIL
The only special thing I have is a maple.ini with kernelopts(floatPi=false):
Changing it to true doesn't make any difference.

@Un4 In Maple 2015 this gives a reasonable plot with 10^6 instead of undefined, which works in Maple 2016:

yield_h := proc (sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u, f__tx, f__ty, alpha)
  local fr:=f__r(sigma__x, sigma__y, tau__xy, f__tx, f__ty, alpha),
        fh:=f__h(sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u);
if evalf(fr) < 0 then fh else 10^6 end if end proc;
yield_r := proc (sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u, f__tx, f__ty, alpha)
  local fr:=f__r(sigma__x, sigma__y, tau__xy, f__tx, f__ty, alpha),
        fh:=f__h(sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u);
if evalf(fh) < 0 then fr else 10^6 end if end proc;
plots:-implicitplot3d([rcurry(yield_h,5$9),rcurry(yield_r,5$9)], -10..10, -10 .. 10, 0 .. 5, style = surfacecontour, grid=[50,50,50], axes = normal,shading=[zhue,zgrayscale);


In Maple 2016 with undefined instead of the number 10^6:

You will notice the difference if you have both versions, turn the graphs and look under the skirt :-)

@Un4 By searching MaplePrimes I found

http://www.mapleprimes.com/posts/137640-Two-Variants

where Kitonum's example

F := piecewise(x >= 0 and y >= 0 and z >= 0, 2*x+3*y+z-6, undefined):
plots:-implicitplot3d(F, x = -3 .. 3, y = -3 .. 3, z = -1 .. 7, color = red, axes = normal, scaling = constrained, style = surface, numpoints = 50000, view=[-3 .. 3, -3 .. 3, -1 .. 7]);

was described by Markiyan Hirnyk as not working in Maple 16.01, but OK in Maple 13.
The example works in Maple 2016, but not in Maple 2015.

@Un4
Note added: You use Maple 2015. I used Maple 2016. When I use Maple 2015 I get an empty plot too with undefined.

I changed as I suggested in my simple example 9999999 to undefined:

yield := proc (sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u, f__tx, f__ty, alpha) if evalf(f__r(sigma__x, sigma__y, tau__xy, f__tx, f__ty, alpha)) < 0 then return f__h(sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u) else undefined end if end proc;

Then tried (all parameters set at 5):
plots:-implicitplot3d('yield(sigma__x, sigma__y, tau__xy,5$9)' = 0, sigma__x = -10 .. 10, sigma__y = -10 .. 10, tau__xy = 0 .. 10, style = surfacecontour, numpoints = 100000, axes = normal);

and I got (in Maple 2016):

@Un4 Could you give those functions as text, not images. Alternatively, upload a worksheet. Typing is not my favorite pastime.

I'm confused about your talk about two surfaces. I only see one on the image.

@Kitonum I never use it either. If you have two algebraic expressions you suspect are equal I find it natural to do as you do: Subtract the two from each other and do some manipulation, in this case expand.

@Kitonum The answer is plain and simple.

In similar but slightly more complicated situations one must be careful:

solve(x^3-3*x+1=0,x); # Each root found by solve has an I.
evalc([%]); #All 3 roots are real.

@Carl Love You mentioned a named expression sequence and gave an example of a list of lists.
So consider:
f:=(A, index::list)-> A[index[]];
LoL:=[[2,3],[8,9]];
map(f,LoL,[-1]); #Works as intended
map(`?[]`, LoL, [-1]); # So does this
SoL:=[2,3],[8,9]; #Named sequence of lists
map(`?[]`, SoL, [-1]); #As expected it doesn't work: Error
map(f, SoL, [-1]); # Doesn't work as maybe intended
## But the following works for `?[]`, but not for f:
`?[]`(SoL,[-1]); #Works
f(SoL,[-1]); #Error

To me this sounds like the familiar shooting method.
Do I understand you correctly if I interpret your statement
   "should work the package of optimization in relation to a system of nonlinear equations"
as meaning that instead of using fsolve (or solve) on the equation y(1) = Y1 you would minimize (e.g.) the (sum of) square(s) y(1) -Y1 ?
That would just amount to using another solver than fsolve, e.g. Optimization:-LSSolve.