// pde44_eq.swift  discretize, build linear equations, solve
// solve  uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
//     4  utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
//     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
//     10 uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =
//     f(x,y,z,t)
//     f(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +
//                  180 z^2*t + 200 y^2z*t + 44 t^3 +
//                  110 y^2z^2t + 66 xy + 12t^4 +
//                  24 z^4 + 36 y^4 + 48 x^4 +
//                  60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t 
//
// boundary conditions computed using ub(x,y,z,t)
// analytic solution is ub(x,y,z,t) = t^4 + 2 z^4 + 3 y^4 + 4 x^4 +
//                                   5 y^2 z^2 t^2 + 6 x y t + 
//                                   7 x z + 8 t + 9
//
// replace continuous derivatives with discrete derivatives
// solve for Uijkl numerical approximation U(x_i,y_j,z_k,t_l)
// A * U = F, solve simultaneous equations for U
// 
//
// fourth order numerical difference order=4, npoints=5, term=0
// d^4u/dx^4 = (1/(b4[0]*hx^4))*(a4[0][0]*u(x,y,z,t) + a4[0][1]*u(x+hx,y,z,t) +
// a4[0][2]*u(x+2hx,y,z,t) + a4[0][3]*u(x+3hx,y,z,t) +a4[0][4]*u(x+4hx,y,z,t))
//
// d^4u/dx^4 using subscripts i,j,k,l for i a minimum subscript of x
// d^4u/dx^4 =(1/(b4[0]*hx^4))*(a4[0][0]*u(i,j,k,l) + a4[0][1]*u(i+1,j,k,l) +
// a4[0][2]*u(i+2,j,k,l) + a4[0][3]*u(i+3,j,k,l) +a4[0][4]*u(i+4,j,k,l))
//
// d^4u/dy^4 using subscripts i,j,k,l for j a minimum subscript of y, etc 
// d^4u/dy^4 =(1/(b4[0]*hy^4))*(a4[0][0]*u(i,j,k,l) + a4[0][1]*u(i,j+1,k,l) +
// a4[0][2]*u(i,j+2,k,l) + a4[0][3]*u(i,j+3,k,l) +a4[0][4]*u(i,j+4,k,l))
//
// d^4u/dx^4 using subscripts i,j,k,l for i-2 and i+2 allowed subscripts of x
// d^4u/dx^4 =(1/(b4[2]*hx^4))*(a4[2][0]*u(i-2,j,k,l) + a4[2][1]*u(i-1,j,k,l) +
// a4[2][2]*u(i,j,k,l) + a4[2][3]*u(i+1,j,k,l) +a4[2][4]*u(i+2,j,k,l))
//
// d^4u/dx^4 using subscripts i,j,k,l for i a maximum subscript of x
// d^4u/dx^4 =(1/(b4[4]*hx^4))*(a4[4][0]*u(i-4,j,k,l) + a4[4][1]*u(i-3,j,k,l) +
// a4[4][2]*u(i-2,j,k,l) + a4[4][3]*u(i-1,j,k,l) +a4[4][4]*u(i,j,k,l))
//
// subscripts in code use x[i], y[ii], z[iii], t[iiii]
// derivatives are computed with b,hx,a included as in
// cx[j] for first derivative of x at point j
// cxxxx[j] for fourth derivative of x at point j
// cyy[jj]*czz[jjj] for d^4u/dy^2dz^2 at jj, jjj
 
// The subscript versions of the derivatives are substituted into the
// equation to be solved. The resulting equation is analytically solved
// for u(i,j,k,l) = other u terms with coeficients and + f(xi,yj,zk,tl)
// The boundaries xmin..xmax, ymin..ymax, zmin..zmax, tmin..tmax are known
// The number of grid points in each dimension is known nx, ny, nz, nt
// hx=(xmax-xmin)/(nx-1)
// hy=(ymax-ymin)/(ny-1)
// hz=(zmax-zmin)/(nz-1)
// ht=(tmax-tmin)/(tx-1) 
// thus xi=xmin+i*hx, yj=ymin+j*hy, zk=zmin+k*hz, tl=tmin+l*ht
// then the linear system of equations is built,
// nxyzt=(nx-2)*(ny-2)*(nz-2)*(nt-2) equations because the boundary value
//                                   is known on each end.
// The matrix A is nxyzt rows by  (nxyzt+1) columns including F
// The forcing function vector F is nxyzt elements adjunct to A
// The solution vector U is nxyzt elements.
// The building of the A matix requires 4 nested loops for rows
// i  in 1..nx-1, j  in 1..ny-1,  k  in 1..nz-1, l  in 1..nt-1 for rows of A
// then each term in the differential equation fills in that row
// and cs = (nx-2)*(ny-2)*(nz-2)*(nt-2) the RHS, fills the last column of A
// subscripting functions are used to compute linear subscripts of A, U
// row=(i-1)*(ny-2)*(nz-2)*(nt-2)+(ii-1)*(nz-2)*(nt-2)+(iii-1)*(nt-2)+(iiii-1)
// row*(nxyzt+1)+col for A
//
// care must be taken to move the negative of boundary elements to
// the right hand size of the equation to be solved. These are subtracted
// from the computed value of the forcing function, cs, for that row.
// tests are j==0 || j==nx-1   ...  jjjj==0 || jjjj==nt-1 for boundaries
 
// uses simeqb should be imported
// uses nuderiv should be imported
// uses inverse should be imported
// also other utility functions for pde44_eq below

print("pde44_eq.swift running")
print("Given uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) + ")
print("    4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) + ")
print("    7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) + ")
print("    10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) = ")
print("  f(x,y,z,t) ")
print("  f(x,y,z,t) = 706 + 120 t^2 + 140 y^2z + 768 x +   ")
print("               180 z^2*t + 200 y^2z*t + 44 t^3 +    ")
print("               110 y^2z^2t + 66 xy + 12t^4 +        ")
print("               24 z^4 + 36 y^4 + 48 x^4 +           ")
print("               60 y^2z^2t^2 + 72 xyt + 84 xz + 96 t ") 
print("xmin<=x<=xmax ymin<=y<=ymax zmin<=z<=zmax Boundaries ")
print("Analytic solution u(x,y,z,t) =  t^4 + 2 z^4 + 3 y^4 + ")
print("     4 x^4 + 5 y^2 z^2 t^2 + 6 x y t + 7 x z + 8 t + 9 ")
print(" ")

let nx = 6 // problem parameters
let ny = 6
let nz = 6
let nt = 6
let nxyzt = (nx-2)*(ny-2)*(nz-2)*(nt-2)
print("nx=\(nx), ny=\(ny), nz=\(nz), nt=\(nt), nxyzt=\(nxyzt)")

var A = [[Double]](repeating:[Double](repeating:0.0,count:nxyzt+1),count:nxyzt)
let cs = (nx-2)*(ny-2)*(nz-2)*(nt-2) // constant RHS column 
                                     // last column is RHS 
var U = [Double](repeating:0.0,count:nxyzt)
var Ua = [Double](repeating:0.0,count:nxyzt)


func s(_ i:Int,_ ii:Int,_ iii:Int,_ iiii:Int) -> Int { // for x,y,z,t
  return (i-1)*(ny-2)*(nz-2)*(nt-2) + (ii-1)*(nz-2)*(nt-2) + (iii-1)*(nt-2) +   (iiii-1)
} // end s

func sk(_ k:Int,_ i:Int,_ ii:Int,_ iii:Int,_ iiii:Int) -> Int { // for x,y,z,t 
  return k*(nxyzt+1) + (i-1)*(ny-2)*(nz-2)*(nt-2) + (ii-1)*(nz-2)*(nt-2) + (iii-1)*(nt-2) + (iiii-1)
} // end sk 

func ss(_ i:Int,_ ii:Int,_ iii:Int,_ iiii:Int,
        _ j:Int,_ jj:Int,_ jjj:Int,_ jjjj:Int) -> Int {
  return s(i,ii,iii,iiii)*(nxyzt+1) + s(j,jj,jjj,jjjj)
} // end ss 

func initialize_A() -> [[Double]] { 
  for i in 1..<nx-1 {
    for ii in 1..<ny-1 {
      for iii in 1..<nz-1 {
        for iiii in 1..<nt-1 {
          for j in 1..<nx-1    {
            for jj in 1..<ny-1   {
              for jjj in 1..<nz-1  {
                for jjjj in 1..<nt-1 {
	          A[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)] = 0.0
	        } // jjjj 
	      } // jjj 
	    } // jj 
          } // j 
        } // iiii 
      } // iii 
    } // ii 
  } // i
  print("A initialized")
  return A
} // end initialize_A

func printA() {
  var k = 0
  for i in 1..<nx-1 {
    for ii in 1..<ny-1 {
      for iii in 1..<nz-1 {
        for iiii in 1..<nt-1 {
          k = s(i,ii,iii,iiii) // row index 
          for j in 1..<nx-1    {
            for jj in 1..<ny-1   {
              for jjj in 1..<nz-1  {
                for jjjj in 1..<nt-1 {
                  print("A[\(i),\(ii),\(iii),\(iiii)][\(j),\(jj),\(jjj),\(jjjj)] = \(A[s(i,ii,iii,iiii)][s(j,jj,jjj,jjjj)])")
	        } // jjjj 
	      } // jjj 
	    } // jj 
          } // j 
          print("RHS A[\(k),\(cs)]=\(A[k][cs])")
          // equivalent A[sk(k,i,ii,iii,iiii+1)] 
        } // iiii 
      } // iii 
    } // ii 
  } // i 
  print(" ")
} // end printA 

// weights of points in each erivative, computed for specific point
var cx = [Double](repeating:0.0,count:nx)
var cy = [Double](repeating:0.0,count:ny)
var cz = [Double](repeating:0.0,count:nx)
var ct = [Double](repeating:0.0,count:nt)
var cxx = [Double](repeating:0.0,count:nx)
var cyy = [Double](repeating:0.0,count:ny)
var czz = [Double](repeating:0.0,count:nz)
var ctt = [Double](repeating:0.0,count:nt)
var cxxx = [Double](repeating:0.0,count:nx)
var cyyy = [Double](repeating:0.0,count:ny)
var czzz = [Double](repeating:0.0,count:nz)
var cttt = [Double](repeating:0.0,count:nt)
var cxxxx = [Double](repeating:0.0,count:nx)
var cyyyy = [Double](repeating:0.0,count:ny)
var czzzz = [Double](repeating:0.0,count:nz)
var ctttt = [Double](repeating:0.0,count:nt) 

let xmin = 0.0 // problem parameters
let xmax = 1.0
let ymin = 0.0
let ymax = 1.0
let zmin = 0.0
let zmax = 1.0
let tmin = 0.0
let tmax = 1.0
print("xmin=\(xmin), xmax=\(xmax), ymin=\(ymin), ymax=\(ymax)")
print("zmin=\(zmin), zmax=\(zmax), tmin=\(tmin), tmax=\(tmax)")

let hx = (xmax-xmin)/(Double)(nx-1)
let hy = (ymax-ymin)/(Double)(ny-1)
let hz = (zmax-zmin)/(Double)(nz-1)
let ht = (tmax-tmin)/(Double)(nt-1)
print("hx=\(hx), hy=\(hy), hz=\(hz), ht=\(ht)")

var coefdxxxx = 0.0 // computed below may not be constant, can use x,y,z,t
var coefdyyyy = 0.0
var coefdzzzz = 0.0
var coefdtttt = 0.0
var coefdxyt = 0.0
var coefdyyzz = 0.0
var coefdztt = 0.0
var coefdxxx = 0.0
var coefdyyt = 0.0
var coefdzt = 0.0
var coefdt = 0.0
var coefu = 0.0

var xg = [Double](repeating:0.0,count:nx)
var yg = [Double](repeating:0.0,count:ny)
var zg = [Double](repeating:0.0,count:nz)
var tg = [Double](repeating:0.0,count:nt)

// initialize grid (does not need to be uniform)
print("x grid and analytic solution at xg[i],ymin,zmin,tmin ")
for i in 0..<nx {
  xg[i] = xmin+(Double)(i)*hx
  Ua[i] = uana(xg[i],ymin,zmin,tmin) // temp 
  print("i=\(i), Ua(\(xg[i]),\(ymin),\(zmin),\(tmin)) = \(Ua[i])")
} // i
print(" ")

print("y grid, and analytic solution at xmin,yg[ii],zmin,tmin ")
for ii in 0..<ny {
  yg[ii] = ymin+(Double)(ii)*hy
  Ua[ii] = uana(xmin,yg[ii],zmin,tmin) // temp 
  print("ii=\(ii), Ua(\(xmin),\(yg[ii]),\(zmin),\(tmin)) = \(Ua[ii])")
} // ii
print(" ")

print("z grid, and analytic solution at xmin,ymin,zg[iii],tmin ")
for iii in 0..<nz {
  zg[iii] = zmin+(Double)(iii)*hz
  Ua[iii] = uana(xmin,ymin,zg[iii],tmin) // temp 
  print("iii=\(iii), Ua(\(xmin),\(ymin),\(zg[iii]),\(tmin)) = \(Ua[iii])")
} // iii
print(" ")

print("t grid, and analytic solution at xmin,ymin,zmin,tg[iiii] ")
for iiii in 0..<nt {
  tg[iiii] = tmin+(Double)(iiii)*ht
  Ua[iiii] = uana(xmin,ymin,zmin,tg[iiii]) // temp 
  print("iii=\(iiii), Ua(\(xmin),\(ymin),\(zmin),\(tg[iiii])) = \(Ua[iiii])")
} // iiii
print(" ")

func eval_derivative(_ xord:Int,_ yord:Int,_ zord:Int,_ tord:Int,_ i:Int,_ ii:Int,_ iii:Int,_ iiii:Int,_ u:[Double]) -> Double {

  let nn = max(nx, max(ny, max(nz, nt)))
  var p1 = [Double](repeating:0.0,count:nn)
  var p2 = [Double](repeating:0.0,count:nn)
  var p3 = [Double](repeating:0.0,count:nn)
  var discrete = 0.0
  let x = xg[i]
  let y = yg[ii]
  let z = zg[iii]
  let t = tg[iiii]

  cx = nuderiv(xord, nx, i, xg)
  cy = nuderiv(yord, ny, ii, yg)
  cz = nuderiv(zord, nz, iii, zg)
  ct = nuderiv(tord, nt, iiii, tg)
  discrete = 0.0

  // four cases of one partial x, y, z, t to any order
  if xord != 0 && yord == 0 && zord == 0 && tord == 0 {
    for j in 0..<nx {
      if j==0 || j==nx-1 { discrete += cx[j]*uana(xg[j],y,z,t) }
      else { discrete += cx[j]*u[s(j,ii,iii,iiii)] }
    }
  }
  else if xord == 0 && yord != 0 && zord == 0 && tord == 0 {
    for jj in 0..<ny {
      if jj==0 || jj==ny-1 { discrete += cy[jj]*uana(x,yg[jj],z,t) }
      else { discrete += cy[jj]*u[s(i,jj,iii,iiii)] }
    }
  }
  else if xord == 0 && yord == 0 && zord != 0 && tord == 0 {
    for jjj in 0..<nz {
      if jjj==0 || jjj==nz-1 { discrete += cz[jjj]*uana(x,y,zg[jjj],t) }
      else { discrete += cz[jjj]*u[s(i,ii,jjj,iiii)] }
    }
  }
  else if xord == 0 && yord == 0 && zord == 0 && tord != 0 {
    for jjjj in 0..<nt {
      if jjjj==0 || jjjj==nt-1 { discrete += ct[jjjj]*uana(x,y,z,tg[jjjj]) }
      else { discrete += ct[jjjj]*u[s(i,ii,iii,jjjj)]  }
    }
  }
  // six cases of two partials xy, xz, xt, yz, yt, zt to any order 
  else if xord != 0 && yord != 0 && zord == 0 && tord == 0 {
    for j in 0..<nx {
      p1[j] = 0.0
      for jj in 0..<ny {
	  if j==0 || j==nx-1 || jj==0 || jj==ny-1 {
	     p1[j] += cy[jj]*uana(xg[j],yg[jj],zg[iii],tg[iiii]) }
          else { p1[j] += cy[jj]*u[s(j,jj,iii,iiii)] }
        discrete += cx[j]*p1[j]
      } // end jj
    } // end j
  }
  else if xord != 0 && yord == 0 && zord != 0 && tord == 0 {
    for j in 0..<nx {
      p1[j] = 0.0
      for jjj in 0..<nz {
	if j==0 || j==nx-1 || jjj==0 || jjj==nz-1 {
	   p1[j] += ct[jjj]*uana(xg[j],yg[ii],zg[jjj],tg[iiii]) }
        else { p1[j] += cz[jjj]*u[s(j,ii,jjj,iiii)] }
        discrete += cx[j]*p1[j]
      } // end jjj
    } // end j
  }
  else if xord != 0 && yord == 0 && zord == 0 && tord != 0 {
    for j in 0..<nx {
      p1[j] = 0.0
      for jjjj in 0..<nt {
        if j==0 || j==nx-1 || jjjj==0 || jjjj==nt-1 {
	   p1[j] += ct[jjjj]*uana(xg[j],yg[ii],zg[iii],tg[jjjj]) }
          else { p1[j] += ct[jjjj]*u[s(j,ii,iii,jjjj)] }
        discrete += cx[j]*p1[j]
      } // end jjjj
    } // end j
  }
  else if xord == 0 && yord != 0 && zord != 0 && tord == 0 {
    for jj in 0..<ny {
      p1[jj] = 0.0
      for jjj in 0..<nz {
        if jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1 {
	   p1[jj] += cz[jjj]*uana(xg[i],yg[jj],zg[jjj],tg[iiii]) }
	  else { p1[jj] += cz[jjj]*u[s(i,jj,jjj,iiii)] }
        discrete += cy[jj]*p1[jj]
      } // end jjj
    } // end jj
  }
  else if xord == 0 && yord != 0 && zord == 0 && tord != 0 {
    for jj in 0..<ny {
      p1[jj] = 0.0
      for jjjj in 0..<nt {
        if jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1 {
           p1[jj] += ct[jjjj]*uana(xg[i],yg[jj],zg[iii],tg[jjjj]) }
        else { p1[jj] += ct[jjjj]*u[s(i,jj,iii,jjjj)] }
        discrete += cy[jj]*p1[jj]
      } // end jjjj
    } // end jj
  }
  else if xord == 0 && yord == 0 && zord != 0 && tord != 0 {
    for jjj in 0..<nz {
      p1[jjj] = 0.0
      for jjjj in 0..<nt {
	if jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1 {
           p1[jjj] += ct[jjjj]*uana(xg[i],yg[ii],zg[jjj],tg[jjjj]) }
        else { p1[jjj] += ct[jjjj]*u[s(i,ii,jjj,jjjj)] }
        discrete += cz[jjj]*p1[jjj]
      } // end jjj
    } // end jjjj
  }
  // four cases of three partials xyz, xyt, xzt, yzt to any order 
  else if xord != 0 && yord != 0 && zord != 0 && tord == 0 {
    for j in 0..<nx {
      p1[j] = 0.0
      for jj in 0..<ny {
        p2[jj] = 0.0
        for jjj in 0..<nz {
	  if j==0 || j==nx-1 || jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1 {
                 p2[jj] += cz[jjj]*uana(xg[j],yg[jj],zg[jjj],tg[iiii]) }
          else { p2[jj] += cz[jjj]*u[s(j,jj,jjj,iiii)] }
          p1[j] += cy[jj]*p2[jj]
	} // end jjj
        discrete += cx[j]*p1[j]
      } // end jj
    } // end j
  }
  else if xord != 0 && yord != 0 && zord == 0 && tord != 0 {
    for j in 0..<nx {
      p1[j] = 0.0
      for jj in 0..<ny {
        p2[jj] = 0.0
        for jjjj in 0..<nt {
	  if j==0 || j==nx-1 || jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1 {
                 p2[jj] += ct[jjjj]*uana(xg[j],yg[jj],zg[iii],tg[jjjj]) }
          else { p2[jj] += ct[jjjj]*u[s(j,jj,iii,jjjj)] }
          p1[j] += cy[jj]*p2[jj]
        } // end jjjj
        discrete += cx[j]*p1[j]
      } // end jj
    } // end j
  }
  else if xord != 0 && yord == 0 && zord != 0 && tord != 0 {
    for j in 0..<nx {
      p1[j] = 0.0
      for jjj in 0..<nz {
        p2[jjj] = 0.0
        for jjjj in 0..<nt {
	  if j==0 || j==nx-1 || jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1 {
                 p2[jjj] += ct[jjjj]*uana(xg[j],yg[ii],zg[jjj],tg[jjjj]) }
          else { p2[jjj] += ct[jjjj]*u[s(j,ii,jjj,jjjj)] }
          p1[j] += cz[jjj]*p2[jjj]
	} // end jjjj
        discrete += cx[j]*p1[j]
      } // end jjj
    } // end j
  }
  else if xord == 0 && yord != 0 && zord != 0 && tord != 0 {
    for jj in 0..<ny {
      p1[jj] = 0.0
      for jjj in 0..<nz {
        p2[jjj] = 0.0
        for jjjj in 0..<nt {
	  if jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1{
             p2[jjj] += ct[jjjj]*uana(xg[i],yg[jj],zg[jjj],tg[jjjj]) }
          else { p2[jjj] += ct[jjjj]*u[s(i,jj,jjj,jjjj)] }
          p1[jj] += cz[jjj]*p2[jjj]
	} // end jjjj
        discrete += cy[jj]*p1[jj]
      } // end jjj
    } // end jj
  }
  // one case of four partials
  else {
    for j in 0..<nx {
      p1[j] = 0.0
      for jj in 0..<ny {
        p2[jj] = 0.0
        for jjj in 0..<nz {
          p3[jjj] = 0.0
          for jjjj in 0..<nt {
            if j==0 || j==nx-1 || jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1 {
              p3[jjj] += ct[jjjj]*uana(xg[j],yg[jj],zg[jjj],tg[jjjj]) }
	    else { p3[jjj] += ct[jjjj]*u[s(j,jj,jjj,jjjj)] }
            p2[jj] += cz[jjj]*p3[jjj]
	  } // end jjjj
          p1[j] += cy[jj]*p2[jj]
	} // end jjj
        discrete += cx[j]*p1[j]
      } // end jj
    } // end j  
  } // end if
  return discrete
} // end eval_deriv  just a test


func check_soln(_ u:[Double]) {
  // solve uxxxx(x,y,z,t) + 2 uyyyy(x,y,z,t) + 3 uzzzz(x,y,z,t) +
  //     4 utttt(x,y,z,t) + 5  uxyt(x,y,z,t) + 6 uyyzz(x,y,z,t) +
  //     7  uztt(x,y,z,t) + 8  uxxx(x,y,z,t) + 9  uyyt(x,y,z,t) +
  //     10  uzt(x,y,z,t) + 11   ut(x,y,z,t) + 12    u(x,y,z,t) =
  //     F(x,y,z,t)
  var val = 0.0
  var err = 0.0
  var smaxerr = 0.0
  for i in 1..<nx-1 {
    for ii in 1..<ny-1 {
      for iii in 1..<nz-1 {
        for iiii in 1..<nt-1 {
          val = 0.0
          // uxxxx(x,y,z,t)
          val += eval_derivative(4, 0, 0, 0, i, ii, iii, iiii, u)
	  // + 2 uyyyy(x,y,z,t)
          val += 2.0*eval_derivative(0, 4, 0, 0, i, ii, iii, iiii, u)
	  // + 3 uzzzz(x,y,z,t)
          val += 3.0*eval_derivative(0, 0, 4, 0, i, ii, iii, iiii, u)
          // + 4 utttt(x,y,z,t)
          val += 4.0*eval_derivative(0, 0, 0, 4, i, ii, iii, iiii, u)
	  // + 5  uxyt(x,y,z,t)
          val += 5.0*eval_derivative(1, 1, 0, 1, i, ii, iii, iiii, u)
	  // + 6 uyyzz(x,y,z,t)
          val += 6.0*eval_derivative(0, 2, 2, 0, i, ii, iii, iiii, u)
          // + 7  uztt(x,y,z,t)
          val += 7.0*eval_derivative(0, 0, 1, 2, i, ii, iii, iiii, u)
	  // + 8  uxxx(x,y,z,t)
          val += 8.0*eval_derivative(3, 0, 0, 0, i, ii, iii, iiii, u)
	  // + 9  uyyt(x,y,z,t)
          val += 9.0*eval_derivative(0, 2, 0, 1, i, ii, iii, iiii, u)
          // +10   uzt(x,y,z,t)
          val += 10.0*eval_derivative(0, 0, 1, 1, i, ii, iii, iiii, u)
	  // +11    ut(x,y,z,t)
          val += 11.0*eval_derivative(0, 0, 0, 1, i, ii, iii, iiii, u)
	  // +12     u(x,y,z,t)
          val += 12.0*uana(xg[i],yg[ii],zg[iii],tg[iiii])
          // -       f(x,y,z,t) should be zero
	  err = abs(val - f(xg[i],yg[ii],zg[iii],tg[iiii]))
          if err>smaxerr {
	    smaxerr = err
	    if err > 1.0 { print("\(i),\(ii),\(iii),\(iiii),err=\(err)") }
	  }
	} // end iiii
      } // end iii
    } // end ii
  } // end i
  print("check_soln maxerr=\(smaxerr)")
} // end check_soln


var val = 0.0
var err = 0.0
var avgerr = 0.0
var maxerr = 0.0

print("check cxxxx = nuderiv(4,nx,0,xg)")
cxxxx = nuderiv(4,nx,0,xg)
for i in 0..<nx { print("cxxxx[\(i)]=\(cxxxx[i])") }
print("check cxxxx = nuderiv(4,nx,nx-1,xg)")
cxxxx = nuderiv(4,nx,nx-1,xg)
for i in 0..<nx { print("cxxxx[\(i)]=\(cxxxx[i])") }
print(" ")

print("put solution in Ua vector ")
// put solution in Ua vector 
var x = 0.0
var y = 0.0
var z = 0.0
var t = 0.0

for i in 1..<nx-1 {
  x = xg[i]
  for ii in 1..<ny-1 {
    y = yg[ii]
    for iii in 1..<nz-1 {
      z = zg[iii]
      for iiii in 1..<nt-1 {
        t = tg[iiii]
        Ua[s(i,ii,iii,iiii)] = uana(x,y,z,t)
      } // iiii 
    } // iii 
  } // ii
} // i 
print(" ")

print("initialize A ")
// initialize A 
A = initialize_A()

var k = 0
// compute entries in A matrix, inside boundary 
print("compute A matrix ")
val = 0.0
for i in 1..<nx-1 {
  for ii in 1..<ny-1 {
    for iii in 1..<nz-1 {
      for iiii in 1..<nt-1 {
        k = s(i,ii,iii,iiii) // row index 

        // for each term in first equation 

	// for d^4u/dx^4 
        coefdxxxx = 1.0
        cxxxx = nuderiv(4,nx,i,xg)
        for j in 0..<nx {
          val = coefdxxxx*cxxxx[j]
          if j==0 || j==nx-1 {
	     A[k][cs] -= val*uana(xg[j],yg[ii],zg[iii],tg[iiii]) }
          else { A[k][s(j,ii,iii,iiii)] += val }
	} // end j

	// for d^4u/dy^4 
        coefdyyyy = 2.0
        cyyyy = nuderiv(4,ny,ii,yg)
        for jj in 0..<ny {
          val = coefdyyyy*cyyyy[jj]
          if jj==0 || jj==ny-1 {
	     A[k][cs] -=  val*uana(xg[i],yg[jj],zg[iii],tg[iiii]) }
          else { A[k][s(i,jj,iii,iiii)] += val }
        } // end jj

	// for d^4u/dz^4 
        coefdzzzz = 3.0
        czzzz = nuderiv(4,nz,iii,zg)
        for jjj in 0..<nz {
          val = coefdzzzz*czzzz[jjj]
          if jjj==0 || jjj==nz-1 {
             A[k][cs] -= val*uana(xg[i],yg[ii],zg[jjj],tg[iiii]) }
          else { A[k][s(i,ii,jjj,iiii)] += val }
        } // end jjj

	// for d^4u/dt^4 
        coefdtttt = 4.0
        ctttt = nuderiv(4,nt,iiii,tg)
        for jjjj in 0..<nt {
          val = coefdtttt*ctttt[jjjj]
          if jjjj==0 || jjjj==nt-1 {
	     A[k][cs] -= val*uana(xg[i],yg[ii],zg[iii],tg[jjjj]) }
          else { A[k][s(i,ii,iii,jjjj)] += val }
	} // end jjjj

	// for d^3u/dxdydt 
        coefdxyt  = 5.0
        cx = nuderiv(1,nx,i,xg)
        cy = nuderiv(1,ny,ii,yg)
        ct = nuderiv(1,nt,iiii,tg)
        for j in 0..<nx {
          for jj in 0..<ny {
            for jjjj in 0..<nt {
              val = coefdxyt*cx[j]*cy[jj]*ct[jjjj]
              if j==0 || j==nx-1 || jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1{
	         A[k][cs] -= val*uana(xg[j],yg[jj],zg[iii],tg[jjjj]) }
              else { A[k][s(j,jj,iii,jjjj)] += val }
	    } // end jjjj 
	  } // end jj 
        } // end j 

	// for d^4u/dy^2dz^2 
        coefdyyzz = 6.0
        cyy = nuderiv(2,ny,ii,yg)
        czz = nuderiv(2,nz,iii,zg)
        for jj in 0..<ny {
          for jjj in 0..<nz {
            val = coefdyyzz*cyy[jj]*czz[jjj]
            if jj==0 || jj==ny-1 || jjj==0 || jjj==nz-1 {
	       A[k][cs] -= val*uana(xg[i],yg[jj],zg[jjj],tg[iiii]) }
            else { A[k][s(i,jj,jjj,iiii)] += val }
	  } // end jjj 
        } // end jj 

	// for d^3u/dzdt^2 
        coefdztt  = 7.0
        cz = nuderiv(1,nz,iii,zg)
        ctt = nuderiv(2,nt,iiii,tg)
        for jjj in 0..<nz {
          for jjjj in 0..<nt {
            val = coefdztt*cz[jjj]*ctt[jjjj]
            if jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1 {
	       A[k][cs] -= val*uana(xg[i],yg[ii],zg[jjj],tg[jjjj]) }
            else { A[k][s(i,ii,jjj,jjjj)] += val }
	  } // end jjjj 
	} // end jjj 

	// for d^3u/dx^3 
        coefdxxx  = 8.0
        cxxx = nuderiv(3,nx,i,xg)
        for j in 0..<nx {
          val = coefdxxx*cxxx[j]
          if j==0 || j==nx-1 {
	     A[k][cs] -= val*uana(xg[j],yg[ii],zg[iii],tg[iiii]) }
          else { A[k][s(j,ii,iii,iiii)] += val }
	} // end j

	// for d^3u/dy^2dt 
        coefdyyt  = 9.0
        cyy = nuderiv(2,ny,ii,yg)
        ct = nuderiv(1,nt,iiii,tg)
        for jj in 0..<ny {
          for jjjj in 0..<nt {
            val = coefdyyt*cyy[jj]*ct[jjjj]
            if jj==0 || jj==ny-1 || jjjj==0 || jjjj==nt-1 {
	       A[k][cs] -= val*uana(xg[i],yg[jj],zg[iii],tg[jjjj]) }
                else { A[k][s(i,jj,iii,jjjj)] += val }
	  } // end jjjj 
	} // end jj 

	// for d^2u/dzdt 
        coefdzt   = 10.0
        cz = nuderiv(1,nz,iii,zg)
        ct = nuderiv(1,nt,iiii,tg)
        for jjj in 0..<nz {
          for jjjj in 0..<nt {
            val = coefdzt*cz[jjj]*ct[jjjj]
            if jjj==0 || jjj==nz-1 || jjjj==0 || jjjj==nt-1 {
	       A[k][cs] -= val*uana(xg[i],yg[ii],zg[jjj],tg[jjjj]) }
            else { A[k][s(i,ii,jjj,jjjj)] += val }
	  } // end jjjj 
	} // end jjj 

	// for du/dt 
        coefdt    = 11.0
        ct = nuderiv(1,nt,iiii,tg)
        for jjjj in 0..<nt {
          val = coefdt*ct[jjjj]
          if jjjj==0 || jjjj==nt-1 {
	     A[k][cs] -= val*uana(xg[i],yg[ii],zg[iii],tg[jjjj]) }
          else { A[k][s(i,ii,iii,jjjj)] += val }
	} // end jjjj

	// for u(x,y,z,t) 
        coefu = 12.0
        A[k][s(i,ii,iii,iiii)] += coefu

        // for f(x,y,z,t) RHS constant 
        A[k][cs] += f(xg[i],yg[ii],zg[iii],tg[iiii])

        // end terms 
      } // end iiii loop 
    } // end iii loop 
  } // end ii loop 
} // end i loop 

// printA()

print("A computed stiffness matrix ")
print(" ")

print("solve \(nxyzt) equations in \(nxyzt) unknowns ")
U = simeqb(A)

print("equations solved ")
print(" ")

print("check PDE against known solution")
check_soln(Ua)
print("check PDE against computed solution")
check_soln(U)
print(" ")

avgerr = 0.0
maxerr = 0.0
print("          U computed, Ua analytic, error ")
for i in 1..<nx-1 {
  for ii in 1..<ny-1 {
    for iii in 1..<nz-1 {
      for iiii in 1..<nt-1 {
        err = abs(U[s(i,ii,iii,iiii)]-Ua[s(i,ii,iii,iiii)])
        if err>maxerr { maxerr = err }
        avgerr = avgerr + err
        print("ug[\(i),\(ii),\(iii),\(iiii)]=\(U[s(i,ii,iii,iiii)]), Ua=\(Ua[s(i,ii,iii,iiii)])), err=\(err)")
      } // iiii 
    } // iii 
  } // ii 
} // i 
avgerr = avgerr/(Double)(nx*ny*nz*nt)
print("                       maxerr=\(maxerr) avgerr=\(avgerr)")
print(" ")
print("end pde44_eq.swift")

func simeqb(_ AY: [[Double]]) -> [Double] { // X 
  let n = AY.count
  let m = n+1
  var X = [Double](repeating:0.0,count:n)
  var B = [[Double]](repeating:[Double](repeating:0.0,count:m),count:n)
  var ROW = [Int](repeating:0,count:n) // ROW INTERCHANGE INDICIES
  var HOLD = 0
  var I_PIVOT = 0      // PIVOT INDICIES
  var PIVOT = 0.0      // PIVOT ELEMENT VALUE
  var ABS_PIVOT = 0.0

  // check sizes
  if m != AY[0].count {
    print("inconistant size arrays rows=\(n),cols=\(AY[0].count)")
    return X
  }
  // BUILD WORKING DATA STRUCTURE 
  for i in 0..<n {
    for j in 0..<m {
      B[i][j] = AY[i][j]
    }
  }

  // SET UP ROW  INTERCHANGE VECTORS
  for k in 0..<n {
    ROW[k] = k
  }

  // BEGIN MAIN REDUCTION LOOP 
  for k in 0..<n {
    // FIND LARGEST ELEMENT FOR PIVOT 
    PIVOT = B[ROW[k]][k]
    ABS_PIVOT = abs(PIVOT)
    I_PIVOT = k
    for i in k..<n {
      if abs(B[ROW[i]][k]) > ABS_PIVOT {
        I_PIVOT = i
        PIVOT = B[ROW[i]][k]
        ABS_PIVOT = abs(PIVOT)
      }
    }

    // HAVE PIVOT, INTERCHANGE ROW POINTERS
    HOLD = ROW[k]
    ROW[k] = ROW[I_PIVOT]
    ROW[I_PIVOT] = HOLD

    // CHECK FOR NEAR SINGULAR
    if ABS_PIVOT < 1.0E-32 {
      for j in k..<m {
        B[ROW[k]][j] = 0.0
      }
      print("redundant row (singular) \(ROW[k])")
    } // singular, delete row
    else {

      // REDUCE ABOUT PIVOT 
      for j in k+1..<m {
        B[ROW[k]][j] = B[ROW[k]][j] / B[ROW[k]][k]
      }

      // INNER REDUCTION LOOP 
      for i in 0..<n {
        if i != k {
          for j in k+1..<m {
            B[ROW[i]][j] = B[ROW[i]][j] - B[ROW[i]][k] * B[ROW[k]][j]
          }
        }
      }
    }
    // FINISHED INNER REDUCTION 
  }

  // END OF MAIN REDUCTION LOOP 
  // BUILD  X  FOR RETURN, UNSCRAMBLING ROWS 
  for i in 0..<n {
    X[i] = B[ROW[i]][n]
  }
  return X
} // end simeqb

// nuderiv should be an import 
// uses  matinv  should be imported from  matrix.swift
func matinv(_ a: [[Double]]) -> [[Double]] {
  let n = a.count
  var inv = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
  for i in 0..<n {
    for j in 0..<n {
      inv[i][j] = a[i][j]
    }
  }

  var row = [Int](repeating:0,count:n)
  var col = [Int](repeating:0,count:n)
  var temp = [Double](repeating:0,count:n)
  var hold = 0
  var I_pivot = 0
  var J_pivot = 0
  var pivot = 0.0
  var abs_pivot = 0.0

  // set up row and column interchange vectors
  for k in 0..<n {
    row[k] = k
    col[k] = k
  }
  // begin main reduction loop
  for k in 0..<n {
    // find largest element for pivot
    pivot = inv[row[k]][col[k]]
    I_pivot = k
    J_pivot = k
    for i in k..<n {
      for j in k..<n {
        abs_pivot = abs(pivot)
        if abs(inv[row[i]][col[j]]) > abs_pivot {
          I_pivot = i
          J_pivot = j
          pivot = inv[row[i]][col[j]]
        }
      }
    }
    if abs(pivot) < 1.0E-10 {
      print("Matrix is singular !")
      return [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
    }
    hold = row[k]
    row[k] = row[I_pivot]
    row[I_pivot] = hold
    hold = col[k]
    col[k] = col[J_pivot]
    col[J_pivot] = hold
    // reduce about pivot
    inv[row[k]][col[k]] = 1.0 / pivot
    for j in 0..<n {
      if j != k {
        inv[row[k]][col[j]] = inv[row[k]][col[j]] * inv[row[k]][col[k]]
      }
    }
    // inner reduction loop
    for i in 0..<n {
      if k != i {
        for j in 0..<n {
          if k != j {
            inv[row[i]][col[j]] = inv[row[i]][col[j]] - inv[row[i]][col[k]] * inv[row[k]][col[j]]
          }
        }
        inv[row[i]][col [k]] = -inv[row[i]][col[k]] * inv[row[k]][col[k]]
      }
    }
  }
  // end main reduction loop

  // unscramble rows
  for j in 0..<n {
    for i in 0..<n {
      temp[col[i]] = inv[row[i]][j]
    }
    for i in 0..<n {
      inv[i][j] = temp[i]
    }
  }
  // unscramble columns
  for i in 0..<n {
    for j in 0..<n {
      temp[row[j]] = inv[i][col[j]]
    }
    for j in 0..<n {
      inv[i][j] = temp[j]
    }
  }
  return inv
} // end matinv   that should be imported

func nuderiv(_ order:Int,_ npoint:Int,_ point:Int,_ x:[Double])->[Double] {
  let n = npoint
  var B = [[Double]](repeating:[Double](repeating:0.0,count:n),count:n)
  var fct = [Double](repeating:0,count:n)
  var c = [Double](repeating:0,count:n)
  var pwr = 0.0

  if order <= 0 { return c }
  for i in 0..<n {   // build matrix that will be inverted 
    pwr = 1.0
    for j in 0..<n {
      B[i][j] = pwr
      pwr = pwr*x[i]
    }
  }
  let A = matinv(B)  // invert the matrix 
  // form derivative for order and point 
  for i in 0..<n { // compute factors for order 
    fct[i] = 1.0
  }
  for j in 0..<order {
    for i in 0..<n {
      fct[i] = fct[i]*(Double)(i-j)
    }
  }
  for i in 0..<n {
    c[i] = 0.0
    pwr = 1.0
    for j in order..<n {
      c[i] = c[i] + fct[j]*A[j][i]*pwr
      pwr = pwr * x[point]
    }
  }
  return c
} // end nuderiv 

// problem boundaries (and analytic solution for checking)
func uana(_ x:Double,_ y:Double,_ z:Double,_ t:Double) -> Double {
  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*y*y*z*z*t*t + 6.0*x*y*t + 7.0*x*z + 8.0*t + 9.0
} // end uana

// problem RHS
func f(_ x:Double,_ y:Double,_ z:Double,_ t:Double) -> Double {
  return 706.0 + 120.0*t*t + 140.0*y*y*z + 768.0*x + 180.0*z*z*t + 200.0*y*y*z*t + 44.0*t*t*t + 110.0*y*y*z*z*t + 66.0*x*y + 12.0*t*t*t*t + 24.0*z*z*z*z + 36.0*y*y*y*y + 48.0*x*x*x*x + 60.0*y*y*z*z*t*t + 72.0*x*y*t + 84.0*x*z + 96.0*t 
} // end f

// end pde44_eq.swift