load("dimension");
  depends(v, [t, x]);
  depends(p, [t, x]);
  /* put known dimensions */
  qput(rho, mass/length^3, dimension);
  qput(v, length/time, dimension);
  qput(t, time, dimension);
  qput(x, length, dimension);
  qput(L, length, dimension);

  /* Navier-Stokes equation */
  ns: rho*(diff(v, t) + v*diff(v, x)) = -diff(p, x) + mu*diff(v, x, 2);

  /* assume dimensions for `p' and `mu' are unknown and derive them */
  eq1: dimension(lhs(ns)) = dimension(part(rhs(ns), 1));
  eq2: dimension(lhs(ns)) = dimension(part(rhs(ns), 2));
  sol: solve([eq1, eq2], [dimension(p), dimension(mu)])[1];

  /* this line should give an error becouse `ns' is dimensionally inconsistent */
  /* dimension(ns); */

  /* put dimensions for `p' and `mu'  */
  put(p, assoc(dimension(p), sol),   'dimension);
  put(mu, assoc(dimension(mu), sol), 'dimension);

  /* now it is OK */
  dimension(ns);

  /* define Reynolds number as dimensionless reverse viscosity */
  one_over_mu_unit: natural_unit(1/mu, [L, rho, v]);
  Re: (1/mu)/one_over_mu_unit;

  print("dimension(mu):          ", dimension(mu));
  print("dimension(p):          ", dimension(p));
  print("Reynolds number:        ", Re);
  print("dimension(ns):        ",   dimension(ns));

mass
dimension(mu):           -----------
                         length time
                            mass
dimension(p):           ------------
                                   2
                        length time
                          rho v L
Reynolds number:         [-------]
                            mu
                           mass
dimension(ns):         -------------
                             2     2
                       length  time