function [ ...
  u, dudt, dudx, dudxx, dudy, dudyy, us, ...
  v, dvdt, dvdx, dvdxx, dvdy, dvdyy, vs, ...
  p, dpdt, dpdx, dpdxx, dpdy, dpdyy, ps ] = all_exppoly ( nu, rho, n, x, y, t )

%*****************************************************************************80
%
%% all_exppoly() evaluates the variables of the exppoly flow.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 August 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Xiaoli Li, Jie Shen,
%    Error analysis of the SAC-MAC scheme for the Navier-Stokes equations,
%    arXiv:1909.05131v1 [math.NA] 8 Sep 2019
%
%  Input:
%
%    real NU, the kinematic viscosity.
%
%    real RHO, the fluid density.
%
%    integer N, the number of nodes.
%
%    real X(N), Y(N), the coordinates of nodes.
%
%    real T, the current time.
%
%  Output:
%
%    real U(N), DUDT(N), DUDX(N), DUDXX(N), DUDY(N), DUDYY(N), US(N),
%    the horizontal velocity values.
%
%    real V(N), DVDT(N), DVDX(N), DVDXX(N), DVDY(N), DVDYY(N), VS(N),
%    the vertical velocity values.
%
%    real P(N), DPDT(N), DPDX(N), DPDXX(N), DPDY(N), DPDYY(N), PS(N),
%    the pressure values.
%
  u =     - exp ( t ) * (        x.^4 -  2.0 * x.^3 +       x.^2 ) .* (  2.0 * y.^3 - 3.0 * y.^2 + y   ) / 256.0;
  dudt =  - exp ( t ) * (        x.^4 -  2.0 * x.^3 +       x.^2 ) .* (  2.0 * y.^3 - 3.0 * y.^2 + y   ) / 256.0;
  dudx =  - exp ( t ) * (  4.0 * x.^3 -  6.0 * x.^2 + 2.0 * x    ) .* (  2.0 * y.^3 - 3.0 * y.^2 + y   ) / 256.0;
  dudxx = - exp ( t ) * ( 12.0 * x.^2 - 12.0 * x    + 2.0        ) .* (  2.0 * y.^3 - 3.0 * y.^2 + y   ) / 256.0;
  dudy =  - exp ( t ) * (        x.^4 -  2.0 * x.^3 +       x.^2 ) .* (  6.0 * y.^2 - 6.0 * y    + 1.0 ) / 256.0;
  dudyy = - exp ( t ) * (        x.^4 -  2.0 * x.^3 +       x.^2 ) .* ( 12.0 * y    - 6.0              ) / 256.0;

  v =       exp ( t ) * (  2.0 * x.^3 - 3.0 * x.^2 + x   ) .* (        y.^4 -  2.0 * y.^3 +       y.^2 ) / 256.0;
  dvdt =    exp ( t ) * (  2.0 * x.^3 - 3.0 * x.^2 + x   ) .* (        y.^4 -  2.0 * y.^3 +       y.^2 ) / 256.0;
  dvdx =    exp ( t ) * (  6.0 * x.^2 - 6.0 * x    + 1.0 ) .* (        y.^4 -  2.0 * y.^3 +       y.^2 ) / 256.0;
  dvdxx =   exp ( t ) * ( 12.0 * x    - 6.0 * x          ) .* (        y.^4 -  2.0 * y.^3 +       y.^2 ) / 256.0;
  dvdy =    exp ( t ) * (  2.0 * x.^3 - 3.0 * x.^2 + x   ) .* (  4.0 * y.^3 -  6.0 * y.^2 + 2.0 * y    ) / 256.0;
  dvdyy =   exp ( t ) * (  2.0 * x.^3 - 3.0 * x.^2 + x   ) .* ( 12.0 * y.^2 - 12.0 * y.^2 + 2.0        ) / 256.0;

  p =       exp ( t )         * ( sin ( pi * y ) - 2.0 / pi );
  dpdt =    exp ( t )         * ( sin ( pi * y ) - 2.0 / pi );
  dpdx =    zeros ( size ( x ) );
  dpdxx =   zeros ( size ( x ) );
  dpdy =    exp ( t ) * pi    *   cos ( pi * y );
  dpdyy = - exp ( t ) * pi.^2 *   sin ( pi * y );

  us = dudt + u .* dudx + v .* dudy + dpdx / rho - nu * ( dudxx + dudyy );
  vs = dvdt + u .* dvdx + v .* dvdy + dpdy / rho - nu * ( dvdxx + dvdyy );
  ps = dudx + dvdy;

  return
end
