restart;        
ContoursWithLabels:= proc(
     Expr::algebraic,
     Range1::(name= range(realcons)), Range2::(name= range(realcons)),
     {contours::{posint, {set,list}(realcons)}:= 8},
     {ImplicitplotOptions::{list,set}({name, name= anything}):= NULL},
     {GraphicOptions::{list,set}({name, name= anything}):= NULL},
     {TextOptions::{list,set}({name, name= anything}):= NULL},
     {Coloring::{list,set}({name, name= anything}):= NULL}
)
local
     r1, r2, f, L1, h, S1, P, r, M, C, T, p, p1, m, n, i,
     x:= lhs(Range1), y:= lhs(Range2)
;
     f:= unapply(Expr, (x,y));
     if contours::posint then
          r1:= rand(convert(rhs(Range1), float));
          r2:= rand(convert(rhs(Range2), float));
          L1:= select(type, (f@op)~({seq([r1,r2](), i= 1..205)}), realcons);
          h:= (L1[-6]-L1[1])/contours;
          S1:= [seq(L1[1]+h/2+h*(n-1), n= 1..contours)]
     else #contours::{set,list}(realcons)
          S1:= [contours[]]
     end if;
     userinfo(1, ContoursWithLabels, print('Contours' = evalf[2](S1)), `\n`);
     r:= k-> rand(20..k-20);
     for C in S1 do
          P:= plots:-implicitplot(
               Expr = C, Range1, Range2,
               gridrefine= 3, ImplicitplotOptions[]
          );
          for p in [plottools:-getdata(P)] do
               p1:= convert(p[3], listlist);
               n:= nops(p1);
               if n < 500 then
                    m:= `if`(40 < n, r(n)(), round(n/2));
                    M[`if`(40 < n, [p1[1..m-11], p1[m+11..n]], [p1])[]]:= NULL;
                    T[[p1[m][], evalf[2](C)]]:= NULL
               else
                    h:= trunc(n/2);
                    m:= r(h)();
                    M[p1[1..m-11], p1[m+11..m+h-11], p1[m+h+11..n]]:= NULL;
                    T[[p1[m][], evalf[2](C)], [p1[m+h][], evalf[2](C)]]:= NULL
               end if
          end do
     end do;
     plots:-display(
          [`if`(
               Coloring = NULL,
               NULL,
               plots:-densityplot(Expr, Range1, Range2, Coloring[])
          ),
          plot([indices(M, 'nolist')], color= black, GraphicOptions[]),
          plots:-textplot([indices(T, 'nolist')], TextOptions[])
         ], 'axes'= 'box', 'gridlines'= false, _rest
     )
end proc:

fcns := {C(eta), T(eta), f(eta)}; bet := 1.0; M := 3.0; pr := 6.0; nb := .2; nt := .7; le := 1.0; bi := 3.0; S := 0.; lam := 1.0; N := 10;
sys := (1+1/bet)*(diff(f(eta), eta, eta, eta))-(M*M)*(diff(f(eta), eta))-2*(diff(f(eta), eta))*(diff(f(eta), eta))+f(eta)*(diff(f(eta), eta, eta)) = 0, diff(T(eta), eta, eta)+pr*(f(eta)*(diff(T(eta), eta))-(diff(T(eta), eta))*T(eta))+pr*nb*(diff(T(eta), eta))*(diff(C(eta), eta))+pr*nt*(diff(T(eta), eta))*(diff(T(eta), eta)) = 0, diff(C(eta), eta, eta)+le*pr*(f(eta)*(diff(C(eta), eta))-(diff(f(eta), eta))*C(eta))+nt*(diff(T(eta), eta, eta))/nb = 0;
bc := f(0) = S, (D(f))(0) = lam, (D(f))(N) = 0, (D(T))(0) = -bi*(1-T(0)), T(N) = 0, C(0) = 1, C(N) = 0;
R := dsolve(eval({bc, sys}), fcns, type = numeric, method = bvp[midrich], output = operator);
psi := unapply(eval(x*f(eta)+T(eta), R), x, eta);

        
ContoursWithLabels(
     psi(eta, x), eta = 0 .. N, x = 0 .. 5, contours=12,
     Coloring= [colorstyle= HUE, colorscheme= ["Red","Yellow","Green","Blue"], style= surface],
     TextOptions= [font= [HELVETICA,BOLD,7], color= black],
     ImplicitplotOptions= [gridrefine= 4],
    
     size= [600,600]
);