Question: How can I add color bar to these 2- and 3-D plot in Maple 2018?

Hi,
I want to have a colorbar aside the curve, but I can't. (My Maple version is 2018).

restart:
with(plots):

# --- constants ---
qc := 0.4:
qh := 0.64:
alpha := 0.8:

# --- parameter ranges ---
deltac_min := 0.01: deltac_max := 0.99:
Tch_min := 0.1:     Tch_max := 1.1:

# --- numerical safety ---
epsDen := 1e-12:
Mcap := 3.0: # cap to avoid infinity

# --- denominator ---
Den := (deltac, Tch) ->
       -2*Tch*deltac*qh - 2*Tch*deltac + 2*Tch*qh
       + 2*deltac*qc + 2*Tch + 2*deltac:

# --- safe functions ---
M1safe := (deltac, Tch) ->
    piecewise(Den(deltac,Tch) > epsDen,
              min(2/sqrt(Den(deltac,Tch)), Mcap),
              Mcap):

M0safe := (deltac, Tch) ->
    piecewise(Den(deltac,Tch) > epsDen,
              min(2*alpha/sqrt(Den(deltac,Tch)), Mcap),
              Mcap):

# --- density plots ---
p1 := densityplot(
       M1safe(deltac,Tch),
       deltac = deltac_min .. deltac_max,
       Tch   = Tch_min .. Tch_max,
       grid  = [200,200],
       style = patchnogrid,
       colorstyle = HUE,
       axes = boxed,
       labels = ["δ_c", "T_ch"],
       title = sprintf("M1(δ_c, T_ch)  (capped at %.2f)", Mcap)
     ):

p2 := densityplot(
       M0safe(deltac,Tch),
       deltac = deltac_min .. deltac_max,
       Tch   = Tch_min .. Tch_max,
       grid  = [200,200],
       style = patchnogrid,
       colorstyle = HUE,
       axes = boxed,
       labels = ["δ_c", "T_ch"],
       title = sprintf("M0(δ_c, T_ch)  (capped at %.2f)", Mcap)
     ):

# --- Contour Den = 0 (to mark physical boundary) ---
cont := contourplot(
          Den(deltac,Tch),
          deltac = deltac_min .. deltac_max,
          Tch   = Tch_min .. Tch_max,
          contours = [0],
          color = black,
          thickness = 2
        ):

# --- Overlay contour ---
p1_final := display(p1, cont):
p2_final := display(p2, cont):

# --- Display side by side ---
display(array([p1_final, p2_final]), scaling = constrained);

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

restart;
with(plots):

# === Constants ===
qc := 0.4:
qh := 0.64:
alpha := 0.8:

# === Denominator ===
Den := (deltac, Tch) ->
     -2*Tch*deltac*qh - 2*Tch*deltac + 2*Tch*qh
     + 2*deltac*qc + 2*Tch + 2*deltac:

# === M1 and M0 ===
M1 := (deltac, Tch) -> piecewise(Den(deltac,Tch)>0, 2/sqrt(Den(deltac,Tch)), undefined):
M0 := (deltac, Tch) -> piecewise(Den(deltac,Tch)>0, 2*alpha/sqrt(Den(deltac,Tch)), undefined):

# === Domain ===
deltac_min := 0.05:  deltac_max := 0.95:
Tch_min := 0.01:     Tch_max := 1.10:   # Physically realistic range

# === 3D Plots ===
p1 := plot3d(M1(deltac,Tch),
             deltac=deltac_min..deltac_max,
             Tch=Tch_min..Tch_max,
             grid=[80,80],
             axes=boxed,
             orientation=[-130,45],
             scaling=constrained,
             style=surfacecontour,
             color=green,
             labels=["δ_c", "T_ch", "M1"],
             labelfont=[TIMES,12],
             title="M1(δ_c, T_ch) for qc=0.4, qh=0.64"):

p2 := plot3d(M0(deltac,Tch),
             deltac=deltac_min..deltac_max,
             Tch=Tch_min..Tch_max,
             grid=[80,80],
             axes=boxed,
             orientation=[-130,45],
             scaling=constrained,
             style=surfacecontour,
             color=blue,
             labels=["δ_c", "T_ch", "M0"],
             labelfont=[TIMES,12],
             title="M0(δ_c, T_ch) for α=0.8, qc=0.4, qh=0.64"):

# === Display vertically (M1 above M0) ===
display(Array(1..2, [p1, p2]), insequence=true, scaling=constrained);
display(Array([[p1], [p2]]));