restart;

#
# The expression gcd(x,y) will return 1,
# for obvious reasons explained in the
# help pages.
#
# Equally obviously, one can 'return' the
# expression unevaluated, by enclosing it
# in unevaluation quotes,
#
# Both are illustrated in the following
#
  a:= gcd(x, y);
  b:= 'gcd(x, y)';
  eval(b);

#
# Now you might be able to control the output
# display by using a command from the
# Typesetting() package. Although my personal
# view is that undocumented commands are often
# undocumented for a reason.
#
# The output might 'look' the way you want, but
# can you ever evaluate/parse/convert it to
# anything meaningful? Does this matter?
#
  c:= Typesetting:-mn(gcd)(x,y);
  eval(c);

Warning, expecting only range variable t in expression exact*solution to be plotted but found names [exact, solution]

restart;

#
# Set up numerical values for all problem parameters
#
  params:=[       k=.5,         tau=.95,      mu=0.1e-1,
                 pi=116.1, vartheta=0.8e-2,  phi=0.25e-2,
            epsilon=0.2e-2,     rho=0.5e-1, beta=0.115e-1,
                chi=0.598e-2,     q=.5,      eta=.2,
              delta=.1,       alpha=0.57e-1,   p=.2,
            Upsilon=1.2
          ]:

#
# Define main function
#
  R:= k*tau*(rho*(Upsilon*(mu+alpha+eta)+chi)+(1-rho)*(Upsilon*(1-q)*eta+mu+beta+chi))
      *
      (pi*(-mu*p+mu+phi)/(mu*(vartheta+mu+phi))+epsilon*pi*(mu*p+vartheta)
      /
      (mu*(vartheta+mu+phi)))/((mu+beta+chi)*(mu+alpha+eta)-chi*eta*(1-q));

#
# Compute "all" derivatives and evaluate numerically.
#
# For the purposes of this calculation "all"
# derivatives, means the derivatives with respect to
# every variable returned by indets(R, name)
#
# Output a list of two element lists where each of
# the latter is
#
# [ varName,
#   eval( diff( R, varName), params )
# ]
#
# [ seq( [j, eval( diff( R, j), params )],j in indets(R, name))];

#
# Compute all "sensitivities" (where the sensitivity
# is as defined in Rouben Rostamian response to the
# OP's earlier post) and evaluate numerically.
#
# For the purposes of this calculation "all" sensitivities
# means the sensitivity with respect to every variable
# returned by indets(R, name)
#
# Output a list of two element lists where each of
# the latter is
#
# [ varName,
#   eval( varName*diff( R, varName)/R, params )
# ]
#
  seq( [j, eval( j*diff( R, j)/R, params )],j in indets(R, name));

#
# Pick one of the variables, define lower
# and upper limits for this variable, along
# with a step size
#
  symb:=beta:
  lsymb:= 0.1:
  usymb:= 0.2:
  tsymb:= 0.01:
#
# Compute the "sensitivity" of all parameters
# for all values of the above selected symbol
#
  params:= [ seq
             ( `if`( lhs(j)=symb,
                     NULL,
                     j
                   ),
              j in params
             )
           ]:
  sens:= [ seq
           ( [ j, [ seq
                    ( [ zz,
                        eval
                        ( j*diff( R, j)/R,
                          [params[], symb=zz]
                        )
                      ],
                      zz = lsymb..usymb, tsymb
                    )
                  ]
             ],
             j in indets(params)
           )
         ]:
#
# Produce all the relevant plots
#
  getPlot:= sym -> seq
                   ( `if`
                     ( j[1]=sym,
                       plots:-pointplot
                              ( j[2],
                                style=line,
                                title=j[1],
                                titlefont=[times, bold, 20],
                                labels=['beta', 'sens'],
                                labelfont=[times, bold, 16],
                                axes=boxed
                              ),
                       NULL
                     ),
                     j in sens
                   ):
  for j in indets(R, name) do
      getPlot(j);
  od;

  restart;
#
# Derive 'simplified' pair of expressions
# from OP's  original input
#
  eqns1:={ 4*x1+7*x2+6*x3 = 186,
          floor((1/2)*x1)+floor((1/5)*x2)+floor((1/3)*x3) = 18,
          floor((1/5)*x1)+floor((1/2)*x2)+floor((1/4)*x3) = 21
         }:
#
# Eliminate x3 from the second and third equations
#
  eqns2:=[ subs( isolate( eqns1[1], x3), eqns1[2]),
           subs( isolate( eqns1[1], x3), eqns1[3])
         ]:
#
# Shift rhs to lhs, to get two expressions
# (which one can test for as zero-valued
#
  eqns3:=[ lhs(eqns2[1])-rhs(eqns2[1]),
           lhs(eqns2[2])-rhs(eqns2[2])
         ]:
#
# Define a function to evaluate the above
# pair of expressions for any supplied arguments
# (We will be looking for the case where both
# expressions return zero
#
  f:=(p,q) -> eval( eqns3, [x1=p, x2=q]):
#
# Specify the limits for x1 and z2
#
  x1lim:=[60, 85]:
  x2lim:=[50, 70]:
#
# For the ranges of values for x1 and x2 defined
# above, evaluate the pair of expressions. If both
# are zero, then retain the argument pair, otherwise
# discard

  ans1:=[ seq
          ( seq
            ( `if`( f(i,j)=[0,0],
                    [i,j],
                    NULL
                  ),
              i = x1lim[1]..x1lim[2], (x1lim[2]-x1lim[1])/501
            ),
            j = x2lim[1]..x2lim[2], (x2lim[2]-x2lim[1])/501
          )
        ]:
#
# Plot all the points for which the solution
# is valid
#
  plots:-pointplot( ans1,
                    view=[ x1lim[1]..x1lim[2],
                           x2lim[1]..x2lim[2]
                         ],
                    labels=[x1, x2],
                    labelfont=[times, bold, 16]
                  );
#
# Use the results derived above and the first
# equation in the original set to get the value
# of x3, for each pair of values x1, x2 in ans1
# and plot the triple in 3d
#
  ans2:=[ seq( [ ans1[j][1],
                 ans1[j][2],
                 rhs
                 ( isolate
                   ( eval
                     ( eqns1[1],
                       [ x1=ans1[j][1],
                         x2=ans1[j][2]
                       ]
                     ),
                     x3
                   )
                 )
               ],
               j=1..numelems(ans1)
             )
        ]:
  plots:-pointplot3d( ans2,
                      view=[ x1lim[1]..x1lim[2],
                             x2lim[1]..x2lim[2],
                             min(seq(ans2[j][3], j=1..numelems(ans2)))..max(seq(ans2[j][3], j=1..numelems(ans2)))
                           ],
                      color=black,
                      labels=[x1, x2, x3],
                      labelfont=[times, bold, 16]
                    );