-- pde45_eq.adb  discretize, build linear equations, solve
--               fourth order PDE in 5 dimensions
-- solve 0.5 Uwwww(w,x,y,z,t) + 1.5 Uwxyz(w,x,y,z,t) + 2.5 Uwwtt(w,x,y,z,t) +
--       Uxxxx(w,x,y,z,t) + 2 Uyyyy(w,x,y,z,t) + 3 Uzzzz(w,x,y,z,t) +
--     4 Utttt(w,x,y,z,t) + 5  Uxyt(w,x,y,z,t) + 6 Uyyzz(w,x,y,z,t) +
--     7  Uztt(w,x,y,z,t) + 8  Uxxx(w,x,y,z,t) + 9  Uyyt(w,x,y,z,t) +
--     10  Uzt(w,x,y,z,t) + 11   Ut(w,x,y,z,t) + 12    U(w,x,y,z,t) =
--     F(w,x,y,z,t)
--     F(w,x,y,z,t) = 216*z^2*t*w^2 + 60*y^2*z^2 + 768*x +
--                    168*y^2*z*w^2 + 35*w*z + 144*w^2*t^2 +
--                    108*t*w + 96*x*z + 77*w*x*y*z + 99*w +
--                    10.5*t + 60*w^4+48*x^4 + 36*y^4 +
--                    24*z^4 + 12*t^4 + 70*w*x*y + 72*y^2*z^2*t^2*w^2 +
--                    132*y^2*z^2*t*w^2 + 240*y^2*z*t*w^2 + 44*t^3 +
--                    84*w*x*y*z*t + 660
--
-- boundary conditions computed using U(w,x,y,z,t)  named  uana below
-- analytic solution is U(w,x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +
--                                     5 w^4 + 6 y^2 z^2 t^2 w^2 +
--                                     7 w x y z t + 8 x z + 9 t w + 10
--

with ada.text_io;
use  ada.text_io;
with ada.numerics;
use  ada.numerics;
with ada.numerics.Long_elementary_Functions;
use  ada.numerics.Long_elementary_Functions;
with real_arrays;
use  real_arrays; -- defines type real
with simeq;
with rderiv;
with calendar;

procedure pde45_Eq is
  nw : integer := 6; -- problem parameters
  nx : integer := 6;
  ny : integer := 6;
  nz : integer := 6;
  nt : integer := 6;
  nwxyzt : integer := (nw-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2);

  A: real_matrix(1..nwxyzt,1..nwxyzt); -- A * U = R
  R : real_vector(1..nwxyzt);  -- RHS known
  U : real_vector(1..nwxyzt);  -- solution being computed
  Ua : real_vector(1..nwxyzt); -- solution for checking
  wg : real_vector(1..nw);
  xg : real_vector(1..nx);
  yg : real_vector(1..ny);
  zg : real_vector(1..nz);
  tg : real_vector(1..nt);
  wmin : real := 0.0; -- problem parameters
  wmax : real := 1.1;
  xmin : real := 0.0;
  xmax : real := 1.2;
  ymin : real := 0.0;
  ymax : real := 1.3;
  zmin : real := 0.0;
  zmax : real := 1.4;
  tmin : real := 0.0;
  tmax : real := 1.5;
  hw, hx, hy, hz, ht : real;
  coefdwwww : real; -- computed below if not constant
  coefdwxyz : real;
  coefdwwtt : real;
  coefdxxxx : real;
  coefdyyyy : real;
  coefdzzzz : real;
  coefdtttt : real;
  coefdxyt  : real;
  coefdyyzz : real;
  coefdztt  : real;
  coefdxxx  : real;
  coefdyyt  : real;
  coefdzt   : real;
  coefdt    : real;
  coefu     : real;

  package real_io is new float_io(real);
  use real_io;

  -- internal functions

  function max(x:integer;y:integer) return integer is
  begin
     if x>y then return x; end if;
     return y;
  end max;

  function s(i:integer; ii:integer; iii:integer; iiii:Integer; iiiii:integer)
            return integer is -- for internal w,x,y,z,t
  begin
     return (i-2)*(nx-2)*(ny-2)*(nz-2)*(nt-2) +
            (ii-2)*(ny-2)*(nz-2)*(nt-2) + (iii-2)*(nz-2)*(nt-2) +
            (iiii-2)*(nt-2) + (iiiii-1);
  end S;

  -- problem RHS
  function F(w:real; x:real; y:real; z:real; t:real) return real is
  begin
    return 216.0*z*z*t*w*w + 60.0*y*y*z*z + 768.0*x +
           168.0*y*y*z*w*w + 35.0*w*z + 144.0*w*w*t*t +
           108.0*t*w + 96.0*x*z + 77.0*w*x*y*z + 99.0*w +
           10.5*t + 60.0*w*w*w*w + 48.0*x*x*x*x + 36.0*y*y*y*y +
           24.0*z*z*z*z + 12.0*t*t*t*t + 70.0*w*x*y + 72.0*y*y*z*z*t*t*w*w +
           132.0*y*y*z*z*t*w*w + 240.0*y*y*z*t*w*w + 44.0*t*t*t +
           84.0*w*x*y*z*t + 660.0;
  end F;

  -- problem boundaries (and analytic solution for checking)
  function uana(w:real; x:real; y:real; z:real; t:real)
           return real is
  begin
        -- t^4 + 2 z^4 + 3 y^4 + 4 x^4 +
        -- 5 w^4 + 6 y^2 z^2 t^2 w^2 +
        -- 7 w x y z t + 8 x z + 9 t w + 10
    return t*t*t*t + 2.0*z*z*z*z + 3.0*y*y*y*y + 4.0*x*x*x*x +
           5.0*w*w*w*w + 6.0*y*y*z*z*t*t*w*w +
           7.0*w*x*y*z*t + 8.0*x*z + 9.0*t*w + 10.0;
  end uana;

  procedure printA is -- optional, much printout
     k : integer; -- no cs, R(k)
  begin
    for i in 2..nw-1 loop
      for ii in 2..nx-1 loop
        for iii in 2..ny-1 loop
          for iiii in 2..nz-1 loop
            for iiiii in 2..nt-1 loop
              k := s(i,ii,iii,iiii,iiiii); -- row index
              for j in 2..nw-1 loop
                for jj in 2..nx-1 loop
                  for jjj in 2..ny-1 loop
                    for jjjj in 2..nz-1 loop
                      for jjjjj in 2..nt-1 loop
                        put_line("A(row("&integer'image(i)&","&
                                 integer'image(ii)&","&integer'image(iii)&
                                 ","&integer'image(iiii)&
                                 ","&integer'image(iiiii)&
                                 "),col("&integer'image(j)&
                                 ","&integer'image(jj)&
                                 ","&integer'image(jjj)&
                                 ","&integer'image(jjjj)&
                                 ","&integer'image(jjjjj)&")) ="&
                                 real'image(A(s(i,ii,iii,iiii,iiiii),
                                              s(j,jj,jjj,jjjj,jjjjj))));
                      end loop; -- jjjjj
                    end loop; -- jjjj
                  end loop; -- jjj
                end loop; -- jj
              end loop; -- j
              put_line("RHS R("&integer'image(k)&")="&real'image(R(k)));
            end loop; -- iiiii
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
    put_line(" ");
  end printA;


  function eval_derivative(word:integer; xord:integer; yord:integer;
                           zord:integer; tord:integer;
                           i:integer; ii:integer; iii:integer; iiii:integer;
                           iiiii:integer; u:real_vector) return real is
    cw : real_vector(1..nw);
    cx : real_vector(1..nx);
    cy : real_vector(1..ny);
    cz : real_vector(1..nz);
    ct : real_vector(1..nt);
    nn : integer := max(nw, max(nx, max(ny, max(nz, nt))));
    p1 : real_vector(1..nn);
    p2 : real_vector(1..nn);
    p3 : real_vector(1..nn);
    discrete : real;
    w, x, y, z, t : real;
    debug : Boolean := false;

  begin
    rderiv(word, nw, i, hw, cw);
    w := wg(i);
    rderiv(xord, nx, ii, hx, cx);
    x := xg(ii);
    rderiv(yord, ny, iii, hy, cy);
    y := yg(iii);
    rderiv(zord, nz, iiii, hz, cz);
    z := zg(iiii);
    rderiv(tord, nt, iiiii, ht, ct);
    t := tg(iiiii);
    discrete := 0.0;
    if s(i,ii,iii,iiii,iiiii)=-1 then debug := true; end if;
    if debug then
      put_line("  word="&integer'image(word)&
               ", xord="&integer'image(xord)&
               ", yord="&integer'image(yord)&
               ", zord="&integer'image(zord)&
               ", tord="&integer'image(tord));
      put_line("      i="&integer'image(i)&
               ",    ii="&integer'image(ii)&
               ",   iii="&integer'image(iii)&
               ",  iiii="&integer'image(iiii)&
               ", iiiii="&integer'image(iiiii));
    end if;

    -- five cases of one partial w, x, y, z, t to any order
    if word/=0 and xord=0 and yord=0 and zord=0 and tord=0 then
      for j in 1..nw loop
        if j=1 or j=nw then
          discrete := discrete + cw(j)*uana(wg(j),x,y,z,t);
        else
          discrete := discrete + cw(j)*u(s(j,ii,iii,iiii,iiiii));
        end if;
      end loop; -- j

    elsif word=0 and xord/=0 and yord=0 and zord=0 and tord=0 then
      for jj in 1..nx loop
        if jj=1 or jj=nx then
          discrete := discrete + cx(jj)*uana(w,xg(jj),y,z,t);
        else
          discrete := discrete + cx(jj)*u(s(i,jj,iii,iiii,iiiii));
        end if;
      end loop; -- jj

    elsif word=0 and xord=0 and yord/=0 and zord=0 and tord=0 then
      for jjj in 1..ny loop
        if jjj=1 or jjj=ny then
          discrete := discrete + cy(jjj)*uana(w,x,yg(jjj),z,t);
        else
          discrete := discrete + cy(jjj)*u(s(i,ii,jjj,iiii,iiiii));
        end if;
      end loop; -- jjj

    elsif word=0 and xord=0 and yord=0 and zord/=0 and tord=0 then
      for jjjj in 1..nz loop
        if jjjj=1 or jjjj=nz then
          discrete := discrete + cz(jjjj)*uana(w,x,y,zg(jjjj),t);
        else
          discrete := discrete + cz(jjjj)*u(s(i,ii,iii,jjjj,iiiii));
        end if;
      end loop; -- jjjj

    elsif word=0 and xord=0 and yord=0 and zord=0 and tord/=0 then
      for jjjjj in 1..nt loop
        if jjjjj=1 or jjjjj=nt then
          discrete := discrete + ct(jjjjj)*uana(w,x,y,z,tg(jjjjj));
        else
          discrete := discrete + ct(jjjjj)*u(s(i,ii,iii,iiii,jjjjj));
        end if;
      end loop; -- jjjjj


    -- ten cases of two partials wx, wy, wz, wt,
    --                           xy, xz, xt, yz, yt, zt to any order
    elsif word/=0 and xord/=0 and yord=0 and zord=0 and tord=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jj in 1..nx loop
          if j=1 or j=nw or jj=1 or jj=nx then
            p1(j) := p1(j) + cx(jj)*uana(wg(j),xg(jj),yg(iii),zg(iiii),tg(iiiii));
          else
            p1(j) := p1(j) + cx(jj)*u(s(j,jj,iii,iiii,iiiii));
          end if;
        end loop; -- jj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord=0 and yord/=0 and zord=0 and tord=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jjj in 1..ny loop
          if j=1 or j=nw or jjj=1 or jjj=ny then
            p1(j) := p1(j) + cy(jjj)*uana(wg(j),xg(ii),yg(jjj),zg(iiii),tg(iiiii));
          else
            p1(j) := p1(j) + cy(jjj)*u(s(j,ii,jjj,iiii,iiiii));
          end if;
        end loop; -- jjj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord=0 and yord=0 and zord/=0 and tord=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jjjj in 1..nz loop
          if j=1 or j=nw or jjjj=1 or jjjj=nz then
            p1(j) := p1(j) + cz(jjjj)*uana(wg(j),xg(ii),yg(iii),zg(jjjj),tg(iiiii));
          else
            p1(j) := p1(j) + cz(jjjj)*u(s(j,ii,iii,jjjj,iiiii));
          end if;
        end loop; -- jjjj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord=0 and yord=0 and zord=0 and tord/=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jjjjj in 1..nt loop
          if j=1 or j=nw or jjjjj=1 or jjjjj=nt then
            p1(j) := p1(j) + ct(jjjjj)*uana(wg(j),xg(ii),yg(iii),zg(iiii),tg(jjjjj));
          else
            p1(j) := p1(j) + ct(jjjjj)*u(s(j,ii,iii,iiii,jjjjj));
          end if;
        end loop; -- jjjjj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word=0 and xord/=0 and yord/=0 and zord=0 and tord=0 then
      for jj in 1..nx loop
          p1(jj) := 0.0;
          for jjj in 1..ny loop
          if jj=1 or jj=nx or jjj=1 or jjj=ny then
            p1(jj) := p1(jj) + cy(jjj)*uana(wg(i),xg(jj),yg(jjj),zg(iiii),tg(iiiii));
          else
            p1(jj) := p1(jj) + cy(jjj)*u(s(i,jj,jjj,iiii,iiiii));
          end if;
          end loop; -- jjj
          discrete := discrete + cx(jj)*p1(jj);
       end loop; -- jj


    elsif word=0 and xord/=0 and yord=0 and zord/=0 and tord=0 then
      for jj in 1..nx loop
        p1(jj) := 0.0;
        for jjjj in 1..nz loop
          if jj=1 or jj=nx or jjjj=1 or jjjj=nz then
            p1(jj) := p1(jj) + cz(jjjj)*uana(wg(i),xg(jj),yg(iii),zg(jjjj),tg(iiiii));
          else
            p1(jj) := p1(jj) + cz(jjjj)*u(s(i,jj,iii,jjjj,iiiii));
          end if;
        end loop; -- jjjj
        discrete := discrete + cx(jj)*p1(jj);
      end loop; -- jj

    elsif word=0 and xord/=0 and yord=0 and zord=0 and tord/=0 then
      for jj in 1..nx loop
        p1(jj) := 0.0;
        for jjjjj in 1..nt loop
          if jj=1 or jj=nx or jjjjj=1 or jjjjj=nt then
            p1(jj) := p1(jj) + ct(jjjjj)*uana(wg(i),xg(jj),yg(iii),zg(iiii),tg(jjjjj));
          else
            p1(jj) := p1(jj) + ct(jjjjj)*u(s(i,jj,iii,iiii,jjjjj));
          end if;
        end loop; -- jjjjj
        discrete := discrete + cx(jj)*p1(jj);
      end loop; -- jj

    elsif word=0 and xord=0 and yord/=0 and zord/=0 and tord=0 then
      for jjj in 1..ny loop
        p1(jjj) := 0.0;
        for jjjj in 1..nz loop
          if jjj=1 or jjj=ny or jjjj=1 or jjjj=nz then
            p1(jjj) := p1(jjj) + cz(jjjj)*uana(wg(i),xg(ii),yg(jjj),zg(jjjj),tg(iiiii));
          else
            p1(jjj) := p1(jjj) + cz(jjjj)*u(s(i,ii,jjj,jjjj,iiiii));
          end if;
        end loop; -- jjjj
        discrete := discrete + cy(jjj)*p1(jjj);
      end loop; -- jjj

    elsif word=0 and xord=0 and yord/=0 and zord=0 and tord/=0 then
      for jjj in 1..ny loop
        p1(jjj) := 0.0;
        for jjjjj in 1..nt loop
          if jjj=1 or jjj=ny or jjjjj=1 or jjjjj=nt then
            p1(jjj) := p1(jjj) + ct(jjjjj)*uana(wg(i),xg(ii),yg(jjj),zg(iiii),tg(jjjjj));
          else
            p1(jjj) := p1(jjj) + ct(jjjjj)*u(s(i,ii,jjj,iiii,jjjjj));
          end if;
        end loop; -- jjjjj
        discrete := discrete + cy(jjj)*p1(jjj);
      end loop; -- jjj

    elsif word=0 and xord=0 and yord=0 and zord/=0 and tord/=0 then
      for jjjj in 1..nz loop
        p1(jjjj) := 0.0;
        for jjjjj in 1..nt loop
          if jjjj=1 or jjjj=nz or jjjjj=1 or jjjjj=nt then
            p1(jjjj) := p1(jjjj) + ct(jjjjj)*uana(wg(i),xg(ii),yg(iii),zg(jjjj),tg(jjjjj));
          else
            p1(jjjj) := p1(jjjj) + ct(jjjjj)*u(s(i,ii,iii,jjjj,jjjjj));
          end if;
        end loop; -- jjjjj
        discrete := discrete + cz(jjjj)*p1(jjjj);
      end loop; -- jjjj



    -- ten cases of three partials wxy, wxz, wxt, wyz, wyt, wzt,
    --                             xyz, xyt, xzt, yzt to any order
    elsif word/=0 and xord/=0 and yord/=0 and zord=0 and tord=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jj in 1..nx loop
          p2(jj) := 0.0;
          for jjj in 1..ny loop
            if j=1 or j=nw or
               jj=1 or jj=nx or
               jjj=1 or jjj=ny then
              p2(jj) := p2(jj) + cy(jjj)*uana(wg(j),xg(jj),yg(jjj),zg(iiii),tg(iiiii));
            else
              p2(jj) := p2(jj) + cy(jjj)*u(s(j,jj,jjj,iiii,iiiii));
            end if;
          end loop; -- jjj
          p1(j) := p1(j) + cx(jj)*p2(jj);
        end loop; -- jj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord/=0 and yord=0 and zord/=0 and tord=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jj in 1..nx loop
          p2(jj) := 0.0;
          for jjjj in 1..nz loop
            if j=1 or j=nw or
               jj=1 or jj=nx or
               jjjj=1 or jjjj=nz then
              p2(jj) := p2(jj) + cz(jjjj)*uana(wg(j),xg(jj),yg(iii),zg(jjjj),tg(iiiii));
            else
              p2(jj) := p2(jj) + cz(jjjj)*u(s(j,jj,iii,jjjj,iiiii));
            end if;
          end loop; -- jjjj
          p1(j) := p1(j) + cx(jj)*p2(jj);
        end loop; -- jj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord/=0 and yord=0 and zord=0 and tord/=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jj in 1..nx loop
          p2(jj) := 0.0;
          for jjjjj in 1..nt loop
            if j=1 or j=nw or
               jj=1 or jj=nx or
               jjjjj=1 or jjjjj=nt then
              p2(jj) := p2(jj) + ct(jjjjj)*uana(wg(j),xg(jj),yg(iii),zg(iiii),tg(jjjjj));
            else
              p2(jj) := p2(jj) + ct(jjjjj)*u(s(j,jj,iii,iiii,jjjjj));
            end if;
          end loop; -- jjjjj
          p1(j) := p1(j) + cx(jj)*p2(jj);
        end loop; -- jj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord=0 and yord/=0 and zord/=0 and tord=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jjj in 1..ny loop
          p2(jjj) := 0.0;
          for jjjj in 1..nz loop
            if j=1 or j=nw or
               jjj=1 or jjj=ny or
               jjjj=1 or jjjj=nz then
              p2(jjj) := p2(jjj) + cz(jjjj)*uana(wg(j),xg(ii),yg(jjj),zg(jjjj),tg(iiiii));
            else
              p2(jjj) := p2(jjj) + cz(jjjj)*u(s(j,ii,jjj,jjjj,iiiii));
            end if;
          end loop; -- jjjj
          p1(j) := p1(j) + cy(jjj)*p2(jjj);
        end loop; -- jjj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord=0 and yord/=0 and zord=0 and tord/=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jjj in 1..ny loop
          p2(jjj) := 0.0;
          for jjjjj in 1..nt loop
            if j=1 or j=nw or
               jjj=1 or jjj=ny or
               jjjjj=1 or jjjjj=nt then
              p2(jjj) := p2(jjj) + ct(jjjjj)*uana(wg(j),xg(ii),yg(jjj),zg(iiii),tg(jjjjj));
            else
              p2(jjj) := p2(jjj) + ct(jjjjj)*u(s(j,ii,jjj,iiii,jjjjj));
            end if;
          end loop; -- jjjjj
          p1(j) := p1(j) + cy(jjj)*p2(jjj);
        end loop; -- jjj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord=0 and yord=0 and zord/=0 and tord/=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jjjj in 1..nz loop
          p2(jjjj) := 0.0;
          for jjjjj in 1..nt loop
            if j=1 or j=nw or
               jjjj=1 or jjjj=nz or
               jjjjj=1 or jjjjj=nt then
              p2(jjjj) := p2(jjjj) + ct(jjjjj)*uana(wg(j),xg(ii),yg(iii),zg(jjjj),tg(jjjjj));
            else
              p2(jjjj) := p2(jjjj) + ct(jjjjj)*u(s(j,ii,iii,jjjj,jjjjj));
            end if;
          end loop; -- jjjjj
          p1(j) := p1(j) + cz(jjjj)*p2(jjjj);
        end loop; -- jjjj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word=0 and xord/=0 and yord/=0 and zord/=0 and tord=0 then
      for jj in 1..nx loop
        p1(jj) := 0.0;
        for jjj in 1..ny loop
          p2(jjj) := 0.0;
          for jjjj in 1..nz loop
            if jj=1 or jj=nx or
               jjj=1 or jjj=ny or
               jjjj=1 or jjjj=nz then
              p2(jjj) := p2(jjj) + cz(jjjj)*uana(wg(i),xg(jj),yg(jjj),zg(jjjj),tg(iiiii));
            else
              p2(jjj) := p2(jjj) + cz(jjjj)*u(s(i,jj,jjj,jjjj,iiiii));
            end if;
          end loop; -- jjjj
          p1(jj) := p1(jj) + cy(jjj)*p2(jjj);
        end loop; -- jjj
        discrete := discrete + cx(jj)*p1(jj);
      end loop; -- jj

    elsif word=0 and xord/=0 and yord/=0 and zord=0 and tord/=0 then
      for jj in 1..nx loop
        p1(jj) := 0.0;
        for jjj in 1..ny loop
          p2(jjj) := 0.0;
          for jjjjj in 1..nt loop
            if jj=1 or jj=nx or
               jjj=1 or jjj=ny or
               jjjjj=1 or jjjjj=nt then
              p2(jjj) := p2(jjj) + ct(jjjjj)*uana(wg(i),xg(jj),yg(jjj),zg(iiii),tg(jjjjj));
            else
              p2(jjj) := p2(jjj) + ct(jjjjj)*u(s(i,jj,jjj,iiii,jjjjj));
            end if;
          end loop; -- jjjjj
          p1(jj) := p1(jj) + cy(jjj)*p2(jjj);
        end loop; -- jjj
        discrete := discrete + cx(jj)*p1(jj);
      end loop; -- jj

    elsif word=0 and xord/=0 and yord=0 and zord/=0 and tord/=0 then
      for jj in 1..nx loop
        p1(jj) := 0.0;
        for jjjj in 1..nz loop
          p2(jjjj) := 0.0;
          for jjjjj in 1..nt loop
            if jj=1 or jj=nx or
               jjjj=1 or jjjj=nz or
               jjjjj=1 or jjjjj=nt then
              p2(jjjj) := p2(jjjj) + ct(jjjjj)*uana(wg(i),xg(jj),yg(iii),zg(jjjj),tg(jjjjj));
            else
              p2(jjjj) := p2(jjjj) + ct(jjjjj)*u(s(i,jj,iii,jjjj,jjjjj));
            end if;
          end loop; -- jjjjj
          p1(jj) := p1(jj) + cz(jjjj)*p2(jjjj);
        end loop; -- jjjj
        discrete := discrete + cx(jj)*p1(jj);
      end loop; -- jj

    elsif word=0 and xord=0 and yord/=0 and zord/=0 and tord/=0 then
      for jjj in 1..ny loop
        p1(jjj) := 0.0;
        for jjjj in 1..nz loop
          p2(jjjj) := 0.0;
          for jjjjj in 1..nt loop
            if jjj=1 or jjj=ny or
               jjjj=1 or jjjj=nz or
               jjjjj=1 or jjjjj=nt then
              p2(jjjj) := p2(jjjj) + ct(jjjjj)*uana(wg(i),xg(ii),yg(jjj),zg(jjjj),tg(jjjjj));
            else
              p2(jjjj) := p2(jjjj) + ct(jjjjj)*u(s(i,ii,jjj,jjjj,jjjjj));
            end if;
          end loop; -- jjjjj
          p1(jjj) := p1(jjj) + cz(jjjj)*p2(jjjj);
        end loop; -- jjjj
        discrete := discrete + cy(jjj)*p1(jjj);
      end loop; -- jjj

    -- five cases of four partials wxyz, wxyt, wxzt, wyzt, xyzt
    elsif word/=0 and xord/=0 and zord/=0 and zord/=0 and tord=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jj in 1..nx loop
          p2(jj) := 0.0;
          for jjj in 1..ny loop
            p3(jjj) := 0.0;
            for jjjj in 1..nz loop
              if   j=1 or j=nw or
                   jj=1 or jj=nx or
                   jjj=1 or jjj=ny or
                   jjjj=1 or jjjj=nz then
                p3(jjj) := p3(jjj) + cz(jjjj)*uana(wg(j),xg(jj),yg(jjj),zg(jjjj),tg(iiiii));
              else
                p3(jjj) :=  p3(jjj) + cz(jjjj)*u(s(j,jj,jjj,jjjj,iiiii));
              end if;
            end loop; -- jjjj
            p2(jj) := p2(jj) + cy(jjj)*p3(jjj);
          end loop; -- jjj
          p1(j) := p1(j) + cx(jj)*p2(jj);
        end loop; -- jj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord/=0 and yord/=0 and zord=0 and tord/=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jj in 1..nx loop
          p2(jj) := 0.0;
          for jjj in 1..ny loop
            p3(jjj) := 0.0;
            for jjjjj in 1..nt loop
              if   j=1 or j=nw or
                   jj=1 or jj=nx or
                   jjj=1 or jjj=ny or
                   jjjjj=1 or jjjjj=nt then
                p3(jjj) := p3(jjj) + ct(jjjjj)*uana(wg(j),xg(jj),yg(jjj),zg(iiii),tg(jjjjj));
              else
                p3(jjj) :=  p3(jjj) + ct(jjjjj)*u(s(j,jj,jjj,iiii,jjjjj));
              end if;
            end loop; -- jjjjj
            p2(jj) := p2(jj) + cy(jjj)*p3(jjj);
          end loop; -- jjj
          p1(j) := p1(j) + cx(jj)*p2(jj);
        end loop; -- jj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord/=0 and yord=0 and zord/=0 and tord/=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jj in 1..nx loop
          p2(jj) := 0.0;
          for jjjj in 1..nz loop
            p3(jjjj) := 0.0;
            for jjjjj in 1..nt loop
              if   j=1 or j=nw or
                   jj=1 or jj=nx or
                   jjjj=1 or jjjj=nz or
                   jjjjj=1 or jjjjj=nt then
                p3(jjjj) := p3(jjjj) + ct(jjjjj)*uana(wg(j),xg(jj),yg(iii),zg(jjjj),tg(jjjjj));
              else
                p3(jjjj) :=  p3(jjjj) + ct(jjjjj)*u(s(j,jj,iii,jjjj,jjjjj));
              end if;
            end loop; -- jjjjj
            p2(jj) := p2(jj) + cz(jjjj)*p3(jjjj);
          end loop; -- jjjj
          p1(j) := p1(j) + cx(jj)*p2(jj);
        end loop; -- jj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word/=0 and xord=0 and yord/=0 and zord/=0 and tord/=0 then
      for j in 1..nw loop
        p1(j) := 0.0;
        for jjj in 1..ny loop
          p2(jjj) := 0.0;
          for jjjj in 1..nz loop
            p3(jjjj) := 0.0;
            for jjjjj in 1..nt loop
              if   j=1 or j=nw or
                   jjj=1 or jjj=ny or
                   jjjj=1 or jjjj=nz or
                   jjjjj=1 or jjjjj=nt then
                p3(jjjj) := p3(jjjj) + ct(jjjjj)*uana(wg(j),xg(ii),yg(jjj),zg(jjjj),tg(jjjjj));
              else
                p3(jjjj) :=  p3(jjjj) + ct(jjjjj)*u(s(j,ii,jjj,jjjj,jjjjj));
              end if;
            end loop; -- jjjjj
            p2(jjj) := p2(jjj) + cz(jjjj)*p3(jjjj);
          end loop; -- jjjj
          p1(j) := p1(j) + cy(jjj)*p2(jjj);
        end loop; -- jjj
        discrete := discrete + cw(j)*p1(j);
      end loop; -- j

    elsif word=0 and xord/=0 and yord/=0 and zord/=0 and tord/=0 then
      for jj in 1..nx loop
        p1(jj) := 0.0;
        for jjj in 1..ny loop
          p2(jjj) := 0.0;
          for jjjj in 1..nz loop
            p3(jjjj) := 0.0;
            for jjjjj in 1..nt loop
              if   jj=1 or jj=nx or
                   jjj=1 or jjj=ny or
                   jjjj=1 or jjjj=nz or
                   jjjjj=1 or jjjjj=nt then
                p3(jjjj) := p3(jjjj) + ct(jjjjj)*uana(wg(i),xg(jj),yg(jjj),zg(jjjj),tg(jjjjj));
              else
                p3(jjjj) :=  p3(jjjj) + ct(jjjjj)*u(s(i,jj,jjj,jjjj,jjjjj));
              end if;
            end loop; -- jjjjj
            p2(jjj) := p2(jjj) + cz(jjjj)*p3(jjjj);
          end loop; -- jjjj
          p1(jj) := p1(jj) + cy(jjj)*p2(jjj);
        end loop; -- jjj
        discrete := discrete + cx(jj)*p1(jj);
      end loop; -- jj

    end if; -- on huge list of  elsif
    if debug then
      put_line("discrete derivative = "&real'image(discrete));
    end if;
    return discrete;
  end eval_derivative;

  procedure check_soln(U:real_vector) is
    -- solve 0.5Uwwww(w,x,y,z,t) + 1.5Uwxyz(w,x,y,z,t) + 2.5Uwwtt(w,x,y,z,t) +
    --          Uxxxx(w,x,y,z,t) + 2  Uyyyy(w,x,y,z,t) + 3  Uzzzz(w,x,y,z,t) +
    --        4 Utttt(w,x,y,z,t) + 5   Uxyt(w,x,y,z,t) + 6  Uyyzz(w,x,y,z,t) +
    --        7  Uztt(w,x,y,z,t) + 8   Uxxx(w,x,y,z,t) + 9   Uyyt(w,x,y,z,t) +
    --        10  Uzt(w,x,y,z,t) + 11    Ut(w,x,y,z,t) + 12     U(w,x,y,z,t) =
    --     F(w,x,y,z,t)
    val, err, smaxerr, savgerr : real;
    ierr : integer;
  begin
    smaxerr := 0.0;
    savgerr := 0.0;
    ierr := 0;

    for i in 2..nw-1 loop
      for ii in 2..nx-1 loop
        for iii in 2..ny-1 loop
          for iiii in 2..nz-1 loop
            for iiiii in 2..nt-1 loop
      val := 0.0;
      -- 0.5 Uwwww(w,x,y,z,t)
      val := val +  0.5*eval_derivative(4, 0, 0, 0, 0, i,ii,iii,iiii,iiiii,u);
      -- + 1.5 Uwxyz(w,x,y,z,t)
      val := val +  1.5*eval_derivative(1, 1, 1, 1, 0, i,ii,iii,iiii,iiiii,u);
      -- + 2.5 Uwwtt(w,x,y,z,t)
      val := val +  2.5*eval_derivative(2, 0, 0, 0, 2, i,ii,iii,iiii,iiiii,u);
      -- uxxxx(w,x,y,z,t)
      val := val +      eval_derivative(0, 4, 0, 0, 0, i,ii,iii,iiii,iiiii,u);
      -- + 2 Uyyyy(w,x,y,z,t)
      val := val +  2.0*eval_derivative(0, 0, 4, 0, 0, i,ii,iii,iiii,iiiii,u);
      -- + 3 Uzzzz(w,x,y,z,t)
      val := val +  3.0*eval_derivative(0, 0, 0, 4, 0, i,ii,iii,iiii,iiiii,u);
      -- + 4 Utttt(w,x,y,z,t)
      val := val +  4.0*eval_derivative(0, 0, 0, 0, 4, i,ii,iii,iiii,iiiii,u);
      -- + 5  Uxyt(w,x,y,z,t)
      val := val +  5.0*eval_derivative(0, 1, 1, 0, 1, i,ii,iii,iiii,iiiii,u);
      -- + 6 Uyyzz(w,x,y,z,t)
      val := val +  6.0*eval_derivative(0, 0, 2, 2, 0, i,ii,iii,iiii,iiiii,u);
      -- + 7  Uztt(w,x,y,z,t)
      val := val +  7.0*eval_derivative(0, 0, 0, 1, 2, i,ii,iii,iiii,iiiii,u);
      -- + 8  Uxxx(w,x,y,z,t)
      val := val +  8.0*eval_derivative(0, 3, 0, 0, 0, i,ii,iii,iiii,iiiii,u);
      -- + 9  Uyyt(w,x,y,z,t)
      val := val +  9.0*eval_derivative(0, 0, 2, 0, 1, i,ii,iii,iiii,iiiii,u);
      -- +10   Uzt(w,x,y,z,t)
      val := val + 10.0*eval_derivative(0, 0, 0, 1, 1, i,ii,iii,iiii,iiiii,u);
      -- +11    Ut(w,x,y,z,t)
      val := val + 11.0*eval_derivative(0, 0, 0, 0, 1, i,ii,iii,iiii,iiiii,u);
      -- +12     U(w,x,y,z,t)
      val := val + 12.0*uana(wg(i),xg(ii),yg(iii),zg(iiii),tg(iiiii));
      -- -       F(w,x,y,z,t) should be zero
      err := abs(val - F(wg(i),xg(ii),yg(iii),zg(iiii),tg(iiiii)));
      savgerr := savgerr + err;
      if err>smaxerr then
         smaxerr := err;
         ierr := s(i,ii,iii,iiii,iiiii);
      end if;
            end loop; --iiiii
          end loop; --iiii
        end loop; --iii
      end loop; -- ii
    end loop; -- i
    savgerr := savgerr/real(s(nw-1, nx-1, ny-1, nz-1, nt-1));
    put_line("check_soln avgerr="&real'image(savgerr)&
             ", at"&integer'image(ierr)&" maxerr="&real'image(smaxerr));
  end check_Soln;

  procedure write_soln(which:integer; tag:string) is
    fout : file_type;
  begin
    create(fout, out_file, "pde45_eq_ada"&tag&".dat");

    put_line("writing pde45_eq_ada"&"tag"&".dat");
    put_line(fout,"5 "&integer'image(nw)&integer'image(nx)&
             integer'image(ny)&integer'image(nz)&integer'image(nt));
    for  i in 1..nw loop
      for ii in 1..nx loop
        for  iii in 1..ny loop
          for iiii in 1..nz loop
            for iiiii in 1..nt loop
              put(fout," ");
              put(fout,wg(i),2,6,0);
              put(fout," ");
              put(fout,xg(ii),2,6,0);
              put(fout," ");
              put(fout,yg(iii),2,6,0);
              put(fout," ");
              put(fout,zg(iiii),2,6,0);
              put(fout," ");
              put(fout,tg(iiiii),2,6,0);
              put(fout," ");
              if which=1 then -- expected solution
                put(fout,uana(wg(i),xg(ii),yg(iii),zg(iiii),tg(iiiii)),4,6,0);
              else -- computed solution
                 if i=1 or i=nw or ii=1 or ii=nx or iii=1 or iii=ny or
                    iiii=1 or iiii=nz or iiiii=1 or iiiii=nt then -- boundary
                  put(fout,uana(wg(i),xg(ii),yg(iii),zg(iiii),tg(iiiii)),4,6,0);
                else
                  put(fout,U(S(i,ii,iii,iiii,iiiii)),4,6,0);
                end if;
              end if;
              new_line(fout);
            end loop; -- iiiii
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
    close(fout);
    put_line("finished writing pde44h_eq_ada"&tag&".dat");
  end write_Soln;

  procedure initializeA is
  begin
    put_line("initialize A ");
    -- initialize A
    for i in 2..nw-1 loop
      for ii in 2..nx-1 loop
        for iii in 2..ny-1 loop
          for iiii in 2..nz-1 loop
            for iiiii in 2..nt-1 loop
              for j in 2..nw-1 loop
                for jj in 2..nx-1 loop
                  for jjj in 2..ny-1 loop
                    for jjjj in 2..nz-1 loop
                      for jjjjj in 2..nt-1 loop
                       A(s(i,ii,iii,iiii,iiiii),s(j,jj,jjj,jjjj,jjjjj)) := 0.0;
                      end loop; -- jjjjj
                    end loop; -- jjjj
                  end loop; -- jjj
                end loop; -- jj
              end loop; -- j
            end loop; -- iiiii
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
  end initializeA;


  procedure buildA is -- system of equations
    val : real := 0.0;
    k: integer;
    cw, cww,       cwwww : real_vector(1..nw);
    cx,      cxxx, cxxxx : real_vector(1..nx);
    cy, cyy,       cyyyy : real_vector(1..ny);
    cz, czz,       czzzz : real_vector(1..nz);
    ct, ctt,       ctttt : real_vector(1..nt);
  begin
    -- compute entries in A matrix, inside boundary
    put_line("compute A matrix and y RHS ");

    -- debug print
    put_line("check rderiv(4,nx,1,hx,cxxxx");
    rderiv(4,nx,1,hx,cxxxx);
    for i in 1..nx loop put_line("cxxxx("&integer'image(i)&
                                 ")="&real'image(cxxxx(i))); end loop;
    put_line("check rderiv(4,nx,nx,hx,cxxxx");
    rderiv(4,nx,nx,hx,cxxxx);
    for i in 1..nx loop put_line("cxxxx("&integer'image(i)&
                                 ")="&real'image(cxxxx(i))); end loop;
    put_line(" ");

    val := 0.0;
    for i in 2..nw-1 loop
      for ii in 2..nx-1 loop
        for iii in 2..ny-1 loop
          for iiii in 2..nz-1 loop
            for iiiii in 2..nt-1 loop
              k := s(i,ii,iii,iiii,iiiii); -- row index
              R(k) := 0.0;

        -- for each term in PDE equation

        -- for d^4u/dw^4
        coefdwwww := 0.5;
        rderiv(4,nw,i,hw,cwwww);
        for j in 1..nw loop
          val := coefdwwww*cwwww(j);
          if j=1 or j=nw then
            R(k) := R(k) - val*uana(wg(j),xg(ii),yg(iii),zg(iiii),tg(iiiii));
          else
            A(k,s(j,ii,iii,iiii,iiiii)) := A(k,s(j,ii,iii,iiii,iiiii)) + val;
          end if;
        end loop; -- jj

        -- for d^4u/dwdxdydz
        coefdwxyz := 1.5;
        rderiv(1,nw,i,hw,cw);
        rderiv(1,nx,ii,hx,cx);
        rderiv(1,ny,iii,hy,cy);
        rderiv(1,nz,iiii,hz,cz);
        for j in 1..nw loop
          for jj in 1..nx loop
            for jjj in 1..ny loop
              for jjjj in 1..nz loop
              val := coefdwxyz*cw(j)*cx(jj)*cy(jjj)*cz(jjjj);
              if j=1 or j=nw or jj=1 or jj=nx or jjj=1 or jjj=ny or
                 jjjj=1 or jjjj=nz then
                R(k) := R(k)-val*uana(wg(j),xg(jj),yg(jjj),zg(jjjj),tg(iiiii));
              else
                A(k,s(j,jj,jjj,jjjj,iiiii)) := A(k,s(j,jj,jjj,jjjj,iiiii))+val;
              end if;
              end loop; -- jjjj
            end loop; -- jjj
          end loop; -- jj
        end loop; -- j

        -- for d^4u/dw^2dt^2
        coefdwwtt := 2.5;
        rderiv(2,nw,i,hw,cww);
        rderiv(2,nt,iiiii,ht,ctt);
        for j in 1..nw loop
          for jjjjj in 1..nt loop
            val := coefdwwtt*cww(j)*ctt(jjjjj);
            if j=1 or j=nw or jjjjj=1 or jjjjj=nt then
              R(k) := R(k) - val*uana(wg(j),xg(ii),yg(iii),zg(iiii),tg(jjjjj));
            else
              A(k,s(j,ii,iii,iiii,jjjjj)) := A(k,s(j,ii,iii,iiii,jjjjj)) + val;
            end if;
          end loop; -- jjjjj
        end loop; -- j

        -- for d^4u/dx^4 (should use eval_derivative)
        coefdxxxx := 1.0;
        rderiv(4,nx,ii,hx,cxxxx);
        for jj in 1..nx loop
          val := coefdxxxx*cxxxx(jj);
          if jj=1 or jj=nx then
            R(k) := R(k) - val*uana(wg(i),xg(jj),yg(iii),zg(iiii),tg(iiiii));
          else
            A(k,s(i,jj,iii,iiii,iiiii)) := A(k,s(i,jj,iii,iiii,iiiii)) + val;
          end if;
        end loop; -- jj

        -- for d^4u/dy^4
        coefdyyyy := 2.0;
        rderiv(4,ny,iii,hy,cyyyy);
        for jjj in 1..ny loop
          val := coefdyyyy*cyyyy(jjj);
          if jjj=1 or jjj=ny then
            R(k) := R(k) - val*uana(wg(i),xg(ii),yg(jjj),zg(iiii),tg(iiiii));
          else
            A(k,s(i,ii,jjj,iiii,iiiii)) := A(k,s(i,ii,jjj,iiii,iiiii)) + val;
          end if;
        end loop; -- jjj

        -- for d^4u/dz^4
        coefdzzzz := 3.0;
        rderiv(4,nz,iiii,hz,czzzz);
        for jjjj in 1..nz loop
          val := coefdzzzz*czzzz(jjjj);
          if jjjj=1 or jjjj=nz then
            R(k) := R(k) - val*uana(wg(i),xg(ii),yg(iii),zg(jjjj),tg(iiiii));
          else
            A(k,s(i,ii,iii,jjjj,iiiii)) := A(k,s(i,ii,iii,jjjj,iiiii)) + val;
          end if;
        end loop; -- jjjj

        -- for d^4u/dt^4
        coefdtttt := 4.0;
        rderiv(4,nt,iiiii,ht,ctttt);
        for jjjjj in 1..nt loop
          val := coefdtttt*ctttt(jjjjj);
          if jjjjj=1 or jjjjj=nt then
            R(k) := R(k) - val*uana(wg(i),xg(ii),yg(iii),zg(iiii),tg(jjjjj));
          else
            A(k,s(i,ii,iii,iiii,jjjjj)) := A(k,s(i,ii,iii,iiii,jjjjj)) + val;
          end if;
        end loop; -- jjjjj

        -- for d^3u/dxdydt
        coefdxyt  := 5.0;
        rderiv(1,nx,ii,hx,cx);
        rderiv(1,ny,iii,hy,cy);
        rderiv(1,nt,iiiii,ht,ct);
        for jj in 1..nx loop
          for jjj in 1..ny loop
            for jjjjj in 1..nt loop
              val := coefdxyt*cx(jj)*cy(jjj)*ct(jjjjj);
              if jj=1 or jj=nx or jjj=1 or jjj=ny or
                jjjjj=1 or jjjjj=nt then
                R(k) := R(k)-val*uana(wg(i),xg(jj),yg(jjj),zg(iiii),tg(jjjjj));
              else
                A(k,s(i,jj,jjj,iiii,jjjjj)) := A(k,s(i,jj,jjj,iiii,jjjjj))+val;
              end if;
            end loop; -- jjjjj
          end loop; -- jjj
        end loop; -- jj

        -- for d^4u/dy^2dz^2
        coefdyyzz := 6.0;
        rderiv(2,ny,iii,hy,cyy);
        rderiv(2,nz,iiii,hz,czz);
        for jjj in 1..ny loop
          for jjjj in 1..nz loop
            val := coefdyyzz*cyy(jjj)*czz(jjjj);
            if jjj=1 or jjj=ny or jjjj=1 or jjjj=nz then
              R(k) := R(k) - val*uana(wg(i),xg(ii),yg(jjj),zg(jjjj),tg(iiiii));
            else
              A(k,s(i,ii,jjj,jjjj,iiiii)) := A(k,s(i,ii,jjj,jjjj,iiiii)) + val;
            end if;
          end loop; -- jjjj
        end loop; -- jjj

        -- for d^3u/dzdt^2
        coefdztt  := 7.0;
        rderiv(1,nz,iiii,hz,cz);
        rderiv(2,nt,iiiii,ht,ctt);
        for jjjj in 1..nz loop
          for jjjjj in 1..nt loop
            val := coefdztt*cz(jjjj)*ctt(jjjjj);
            if jjjj=1 or jjjj=nz or jjjjj=1 or jjjjj=nt then
              R(k) := R(k) - val*uana(wg(i),xg(ii),yg(iii),zg(jjjj),tg(jjjjj));
            else
              A(k,s(i,ii,iii,jjjj,jjjjj)) := A(k,s(i,ii,iii,jjjj,jjjjj)) + val;
            end if;
          end loop; -- jjjjj
        end loop; -- jjjj

        -- for d^3u/dx^3
        coefdxxx  := 8.0;
        rderiv(3,nx,ii,hx,cxxx);
        for jj in 1..nx loop
          val := coefdxxx*cxxx(jj);
          if jj=1 or jj=nx then
            R(k) := R(k) - val*uana(wg(i),xg(jj),yg(iii),zg(iiii),tg(iiiii));
          else
            A(k,s(i,jj,iii,iiii,iiiii)) := A(k,s(i,jj,iii,iiii,iiiii)) + val;
          end if;
        end loop; -- jj

        -- for d^3u/dy^2dt
        coefdyyt  := 9.0;
        rderiv(2,ny,iii,hy,cyy);
        rderiv(1,nt,iiiii,ht,ct);
        for jjj in 1..ny loop
          for jjjjj in 1..nt loop
            val := coefdyyt*cyy(jjj)*ct(jjjjj);
            if jjj=1 or jjj=ny or jjjjj=1 or jjjjj=nt then
              R(k) := R(k) - val*uana(wg(i),xg(ii),yg(jjj),zg(iiii),tg(jjjjj));
            else
              A(k,s(i,ii,jjj,iiii,jjjjj)) := A(k,s(i,ii,jjj,iiii,jjjjj)) + val;
            end if;
          end loop; -- jjjjj
        end loop; -- jjj

        -- for d^2u/dzdt
        coefdzt   := 10.0;
        rderiv(1,nz,iiii,hz,cz);
        rderiv(1,nt,iiiii,ht,ct);
        for jjjj in 1..nz loop
          for jjjjj in 1..nt loop
            val := coefdzt*cz(jjjj)*ct(jjjjj);
            if jjjj=1 or jjjj=nz or jjjjj=1 or jjjjj=nt then
              R(k) := R(k) - val*uana(wg(i),xg(ii),yg(iii),zg(jjjj),tg(jjjjj));
            else
              A(k,s(i,ii,iii,jjjj,jjjjj)) := A(k,s(i,ii,iii,jjjj,jjjjj)) + val;
            end if;
          end loop; -- jjjjj
        end loop; -- jjjj

        -- for du/dt
        coefdt    := 11.0;
        rderiv(1,nt,iiiii,ht,ct);
        for jjjjj in 1..nt loop
          val := coefdt*ct(jjjjj);
          if jjjjj=1 or jjjjj=nt then
            R(k) := R(k) - val*uana(wg(i),xg(ii),yg(iii),zg(iiii),tg(jjjjj));
          else
            A(k,s(i,ii,iii,iiii,jjjjj)) := A(k,s(i,ii,iii,iiii,jjjjj)) + val;
          end if;
        end loop; -- jjjjj

        -- for u(w,x,y,z,t)
        coefu     := 12.0;
        A(k,s(i,ii,iii,iiii,iiiii)) := A(k,s(i,ii,iii,iiii,iiiii)) + coefu;

        -- for f(w,x,y,z,t) RHS constant
        R(k) := R(k) + f(wg(i),xg(ii),yg(iii),zg(iiii),tg(iiiii));

        -- end terms

            end loop; -- iiiii
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
  end BuildA;

  procedure main is
    w, x, y, z, t : real;
    err, avgerr, maxerr : real;
    time_start : duration;
    time_now : duration;
    time_last : duration;

  begin -- main
    put_line("pde45_eq.adb running ");
    put_line("Given PDE: ");
    put_line("0.5 Uwwww(w,x,y,z,t) + 1.5 Uwxyz(w,x,y,z,t) + 2.5 Uwwtt(w,x,y,z,t) +");
    put_line("    Uxxxx(w,x,y,z,t) + 2   Uyyyy(w,x,y,z,t) + 3   Uzzzz(w,x,y,z,t) +");
    put_line("  4 Utttt(w,x,y,z,t) + 5    Uxyt(w,x,y,z,t) + 6   Uyyzz(w,x,y,z,t) +");
    put_line("  7  Uztt(w,x,y,z,t) + 8    Uxxx(w,x,y,z,t) + 9    Uyyt(w,x,y,z,t) +");
    put_line(" 10   Uzt(w,x,y,z,t) + 11     Ut(w,x,y,z,t) + 12      U(w,x,y,z,t) =");
    put_line("F(w,x,y,z,t) ");
    put_line("F(w,x,y,z,t) = 216.0*z*z*t*w*w + 60.0*y*y*z*z + 768.0*x +");
    put_line("               168.0*y*y*z*w*w + 35.0*w*z + 144.0*w*w*t*t +");
    put_line("               108.0*t*w + 96.0*x*z + 77.0*w*x*y*z + 99.0*w +");
    put_line("               10.5*t + 60.0*w*w*w*w + 48.0*x*x*x*x + 36.0*y*y*y*y +");
    put_line("               24.0*z*z*z*z + 12.0*t*t*t*t + 70.0*w*x*y + 72.0*y*y*z*z*t*t*w*w +");
    put_line("               132.0*y*y*z*z*t*w*w + 240.0*y*y*z*t*w*w + 44.0*t*t*t +");
    put_line("               84.0*w*x*y*z*t + 660.0 ");
    put_line(" ");
    put_line("wmin<w<wmax xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries ");
    put_line("Analytic solution U(w,x,y,z,t) = ");
    put_line("         t*t*t*t + 2.0*z*z*z*z + 3.0*y*y*y*y + 4.0*x*x*x*x + ");
    put_line("         5.0*w*w*w*w + 6.0*y*y*z*z*t*t*w*w + ");
    put_line("         7.0*w*x*y*z*t + 8.0*x*z + 9.0*t*w + 10.0 ");

    put_line("wmin="&real'image(wmin)&", wmax="&real'image(wmax));
    put_line("xmin="&real'image(xmin)&", xmax="&real'image(xmax));
    put_line("ymin="&real'image(ymin)&", ymax="&real'image(ymax));
    put_line("zmin="&real'image(zmin)&", zmax="&real'image(zmax));
    put_line("tmin="&real'image(tmin)&", tmax="&real'image(tmax));
    put_line("nw="&Integer'Image(nw)&", nx="&integer'image(nx)&
             ", ny="&integer'image(ny)&", nz="&integer'image(nz)&
             ", nt="&integer'image(nt));
    time_start := calendar.seconds(calendar.clock);
    time_last := time_start;

    hw := (wmax-wmin)/real(nw-1);
    put_line("w grid, hw="&real'image(hw));
    for i in 1..nw loop
      wg(i) := wmin+real(i-1)*hw;
      put_line("wg("&integer'image(i)&")="&real'image(wg(i)));
    end loop; -- i
    put_line(" ");
    hx := (xmax-xmin)/real(nx-1);
    put_line("x grid, hx="&real'image(hx));
    for ii in 1..nx loop
      xg(ii) := xmin+real(ii-1)*hx;
      put_line("xg("&integer'image(ii)&")="&real'image(xg(ii)));
    end loop; -- ii
    put_line(" ");
    hy := (ymax-ymin)/real(ny-1);
    put_line("y grid, hy="&real'image(hy));
    for iii in 1..ny loop
      yg(iii) := ymin+real(iii-1)*hy;
      put_line("yg("&integer'image(iii)&")="&real'image(yg(iii)));
    end loop; -- iii
    put_line(" ");
    hz := (zmax-zmin)/real(nz-1);
    put_line("z grid, hz="&real'image(hz));
    for iiii in 1..nz loop
      zg(iiii) := zmin+real(iiii-1)*hz;
      put_line("zg("&integer'image(iiii)&")="&real'image(zg(iiii)));
    end loop; -- iiii
    put_line(" ");
    ht := (tmax-tmin)/real(nt-1);
    put_line("t grid, ht="&real'image(ht));
    for iiiii in 1..nt loop
      tg(iiiii) := tmin+real(iiiii-1)*ht;
      put_line("tg("&integer'image(iiiii)&")="&real'image(tg(iiiii)));
    end loop; -- iiiii
    put_line(" ");

    put_line("put solution in Ua vector ");
    -- put solution in Ua vector
    for i in 2..nw-1 loop
      w := wg(i);
      for ii in 2..nx-1 loop
        x := xg(ii);
        for iii in 2..ny-1 loop
          y := yg(iii);
          for iiii in 2..nz-1 loop
            z := zg(iiii);
            for iiiii in 2..nz-1 loop
              t := tg(iiiii);
              Ua(s(i,ii,iii,iiii,iiiii)) := uana(w,x,y,z,t);
            end loop; -- iiiii
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i
    time_now := calendar.seconds(calendar.clock);
    put_line("time for previous section = "&Duration'image(time_now-time_last)&
             " seconds ");
    time_last := time_now;
    put_line(" ");

    initializeA;

    time_now := calendar.seconds(calendar.clock);
    put_line("time for previous section = "&Duration'image(time_now-time_last)&
             " seconds ");
    time_last := time_now;
    put_line(" ");

    buildA; -- system of equations

    -- printA; -- optional output

    put_line("A computed stiffness matrix ");
    time_now := calendar.seconds(calendar.clock);
    put_line("time for previous section = "&Duration'image(time_now-time_last)&
             " seconds ");
    time_last := time_now;
    put_line(" ");


    put_line("solve "&integer'image(nwxyzt)&" equations in "&
             integer'image(nwxyzt)&" unknowns ");

    U := simeq(A, R);

    put_line("equations solved ");
    time_now := calendar.seconds(calendar.clock);
    put_line("time for previous section = "&Duration'image(time_now-time_last)&
             " seconds ");
    time_last := time_now;
    put_line(" ");

    put_line("check PDE against known solution");
    check_soln(Ua);
    write_soln(1,"Ua");
    put_line("check PDE against computed solution");
    check_soln(U);
    put_line(" ");
    write_soln(2,"");

    avgerr := 0.0;
    maxerr := 0.0;
    put_line("          U computed,ua analytic, error ");
    for i in 2..nw-1 loop
      for ii in 2..nx-1 loop
        for iii in 2..ny-1 loop
          for iiii in 2..nz-1 loop
            for iiiii in 2..nz-1 loop
            err := abs(U(s(i,ii,iii,iiii,iiiii))-Ua(s(i,ii,iii,iiii,iiiii)));
            if err>maxerr then
              maxerr := err;
            end if;
            avgerr := avgerr + err;
            put("ug("&integer'image(i)&","&integer'image(ii)&
                  ","&integer'image(iii)&","&integer'image(iiii)&
                  ","&integer'image(iiiii)&")=");
            put(U(s(i,ii,iii,iiii,iiiii)),3,6); put(", ua=");
            put(Ua(s(i,ii,iii,iiii,iiiii)),3,6); put(", err=");
            put(err);
            new_line;
            end loop; -- iiiii
          end loop; -- iiii
        end loop; -- iii
      end loop; -- ii
    end loop; -- i

    time_now := calendar.seconds(calendar.clock);
    put_line("total CPU time  = "&Duration'image(time_now-time_start)&
             " seconds");
    put_line(" ");


    put_line("                     maxerr="&real'image(maxerr)&
             ", avgerr="&real'image(avgerr/real(nx*ny*nz*nt)));
    put_line(" ");
  end Main;
begin -- pde45_eq
   main;
   put_line("end pde45_eq.adb");
end pde45_eq;