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

%*****************************************************************************80
%
%% all_cavity evaluates the variables of the cavity flow.
%
%  Discussion:
%
%    This is NOT the standard driven cavity problem, in which the lid
%    boundary condition sets the horizontal velocity to 1, and the
%    vertical velocity to 0.
%
%    f(x) = x^4-2x^3+x^2
%    g(y) = y^4 - y^2
%    F(x) = 0.2x^5 - 0.5x^4 + x^3/3
%    F1(x) = - 4x^6 + 12x^5 - 14x^4 + 8x^3 - 2x^2
%    F2(x) = 0.5 * f(x)^2 = 0.5x^8 +2x^7+3x^6+2x^5+0.5*x^4
%    G1(y) = -24y^5 + 8y^3 -4y
%
%    u =   8 f(x)g'(y)
%    v = - 8 f'(x)g(y)
%    p = 8/RE [F(x)g''(y) + f'(x)g'(y)] + 64 F2(x)*(g(y)g''(y)-g'(y)^2)
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 August 2020
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    A C Kengni Jotsa, V A Pennati,
%    A cost-effective FE method for 2D Navier-Stokes equations,
%    Engineering Applications of Computational Fluid Mechanics,
%    Volume 9, Number 1, pages 66-93, 2015.
%
%    Tien-Mo Shih, C H Tan, B C Hwang,
%    Effects of grid staggering on numerical schemes,
%    International Journal for Numerical Methods of Fluids,
%    Volume 9, Number 2, pages 193-212, February 1989.
%
%  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.  Note that DPDXX and DPDYY are returned as
%    zero because life is short.
%
  fint =  0.2 * x.^5 -  0.5 * x.^4 +       x.^3 / 3.0;
  f =           x.^4 -  2.0 * x.^3 +       x.^2;
  fx =    4.0 * x.^3 -  6.0 * x.^2 + 2.0 * x;
  fxx =  12.0 * x.^2 - 12.0 * x    + 2.0;
  fxxx = 24.0 * x    - 12.0;

  g =           y.^4 -        y.^2;
  gy =    4.0 * y.^3 -  2.0 * y;
  gyy =  12.0 * y.^2 -  2.0;
  gyyy = 24.0 * y;

  u =     8.0 * f .* gy;
  dudt = zeros ( size ( x ) );
  dudx =  8.0 * fx .* gy;
  dudxx = 8.0 * fxx .* gy;
  dudy =  8.0 * f .* gyy;
  dudyy = 8.0 * f .* gyyy;

  v =     - 8.0 * fx .* g;
  dvdt = zeros ( size ( x ) );
  dvdx =  - 8.0 * fxx .* g;
  dvdxx = - 8.0 * fxxx .* g;
  dvdy =  - 8.0 * fx .* gy;
  dvdyy = - 8.0 * fx .* gyy;

  p = ( rho / nu ) * 8.0 * ( ...
    ( 0.2 * x.^5 - 0.5 * x.^4 + x.^3 / 3.0 ) .* ( 12.0 * y.^2 - 2.0 ) ...
    + ( 4.0 * x.^3 - 6.0 * x.^2 + 2.0 * x ) .* ( 4.0 * y.^3 - 2.0 * y ) ) ...
    + 32.0 * ( x.^8 + 4.0 * x.^7 + 6.0 * x.^6 + 4.0 * x.^5 + x.^4 ) ...
    .* ( - 4.0 * y.^6 + 2.0 * y.^4 - 2.0 * y.^2 );

  dpdt = zeros ( size ( x ) );
  dpdx = ( rho / nu ) * 8.0 * ( ...
    ( x.^4 - 2.0 * x.^2 + x.^2 ) .* ( 12.0 * y.^2 - 2.0 ) ...
    + ( 12.0 * x.^2 - 12.0 * x + 2.0 ) .* ( 4.0 * y.^3 - 2.0 * y ) ) ...
    + 32.0 * ( 8.0 * x.^7 + 28.0 * x.^6 + 36.0 * x.^5 + 20.0 * x.^4 + 4.0 * x.^3 ) ...
    .* ( - 4.0 * y.^6 + 2.0 * y.^4 - 2.0 * y.^2 );
  dpdxx = zeros ( size ( x ) );
  dpdy = ( rho / nu ) * 8.0 * ( ...
    ( 0.2 * x.^5 - 0.5 * x.^4 + x.^3 / 3.0 ) .* ( 24.0 * y ) ...
    + ( 4.0 * x.^3 - 6.0 * x.^2 + 2.0 * x ) .* ( 12.0 * y.^2 - 2.0 ) ) ...
    + 32.0 * ( x.^8 + 4.0 * x.^7 + 6.0 * x.^6 + 4.0 * x.^5 + x.^4 ) ...
    .* ( - 24.0 * y.^5 + 8.0 * y.^3 - 4.0 * y );
  dpdyy = zeros ( size ( x ) );

  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
