@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.

@snhaseela2 Solving boundary problems for ode systems can be more challenging than solving initial value problems. So any method that works is worth considering. 

When successful dsolve/numeric/bvp is easier to use.

With Douglas Meade's Shoot package you need to rewrite the system as a first order system and you have to provide a start guess. Finding one that works may not be easy.
With the shooting method I gave you don't have to rewrite the system, but for success you need to supply a start guess to the solver used, in my case fsolve. The method as presented in my answer is not written as a ready made package and may be difficult to understand for a new Maple user. Writing a package is not difficult though. I did, but I'm not convinced that the result is good enough for public consumption.

Even dsolve/numeric/bvp in fact needs an initial "profile". It tries to find one itself. If unsuccessful the user can supply an approximate solution, approxsoln=[f(eta)= ...., theta(eta)= ..., phi(eta)= ... ].

@snhaseela2 You have several problems of syntax:

1. All unknown functions (A,B,C, etc) must appear as A(eta), B(eta), C(eta) etc.

2. You need to use another name than D (e.g. D1) since D is used as the differention operator on functions: Thus D(sin) = cos.
3. You cannot use square brackets, [], instead of parentheses, (), for grouping.
4. You cannot leave out some multiplication sign. In 2D math it is enough with a space, but my advice is to use * always.
As an example you have MB without space or *.

Also clearly beta needs to be given a value.
All these syntax problems must be dealt with. Take a close look at eval(ODE, beta = 1) before proceeding: Does it look right?

Probably what you are using is the shoot procedure made by Douglas Meade, which (reasonably enough) needs a start guess. Try the one I used in my answer - and try my answer too!

@snhaseela2 You may not be aware of odeplot from the plots package. It is very useful.
Here is how it could be used in your case using Tom Leslie's revised version of yours.
Notice that I keep all the 4 results in their entirety (as res[1], res[2] etc. )

for k to 4 do
  res[k]:= dsolve(eval({Eq1, Eq2, Eq3, bcs1}, beta = L[k]),numeric, method = bvp[midrich], maxmesh = 4096)
end do;
#Now plot whatever function you like:
plots:-odeplot(res[1],[[eta,D(f)(eta)],[eta,theta(eta)],[eta,(D@@3)(f)(eta)],[eta,-phi(eta)]],0..5);
#Here is a short animation in beta:
S:=seq(plots:-odeplot(res[k],[[eta,D(f)(eta)],[eta,theta(eta)],[eta,(D@@3)(f)(eta)],[eta,-phi(eta)]],0..5),k=1..4):
plots:-display(S,insequence);

##For information about maxmesh go to ?dsolve,numeric,bvp
##The method used by dsolve is not a shooting method. A shooting method formulates the problem as an initial value problem at some point with some given values and some variable ones. Then it tries to find the correct values of the variable ones. In your case you could try shooting from eta=blt (not eta=0 since f(0)=0 is imposed and f(eta) appears in the denominator for the fourth derivative of f. Try
solve({Eq1, Eq2, Eq3},{diff(f(eta),eta$4),diff(theta(eta),eta,eta),diff(phi(eta),eta,eta)});

You define N twice, but don't give a definition of M.
By the way: Use Vector instead of the deprecated vector.

@iman You are plotting w.r.t. sigma1, but sigma1 is set at 40. Did you mean sigma2?

I didn't try solve for very long, so never got a result from that. Switching to fsolve I tried

S:=proc(s2) local res;
   res:=fsolve(eval({Q1, Q2, Q3, Q4},sigma2=s2));
   if not type(res,set) then print(cat("failure at sigma2=",s2)) else eval(p1^2+q1^2,res) end if
end proc;
Digits:=15:
plot(S,-100..200,adaptive=false,numpoints=100);