Using high order discretization to solve differential equations

The goal is to solve the differential equation

a4(x)*Uxxxx(x) + a3(x)*Uxxx(x) + a2(x)*Uxx(x) + a1(x)*U(x) + a0(x)*U(x) = f(x)

f and a0 through a4 being known, computable, analytic functions of x.
The boundary values are known, the value of U(x) the smallest and
largest value of x for the solution space. The total number of
x points will be denoted nx. The two end points are known, thus 
nx-2 is the degrees of freedom.

The desired solution is to find the values of the unknown function U(x)
at a specific set of nx-2 points along the X axis. The method is to set
up nx-2 equations in nx-2 unknowns using f(x) for the right hand side.
Solving the system of linear equations computes the required
unknown values.

Observe the process for numerically computing derivatives:
Given n discrete values of a function, at n unique points along
the X axis, the weights can be determined that compute the
derivative at each point. The points do not need to be
uniformly spaced. 

See comments in nderiv.c, nderiv.java, nderiv.f90, nderiv.adb, etc.
For non uniform spaced grid, see nuderiv for various languages.

There are n weights per point and n points, thus we place the 
weights in an n by n matrix.

The goal is to write a program to solve the differential equation.
This example uses zero based subscripting.

For example, choose n=7. This will allow exact, within roundoff
error, computation of the derivative at every point,
for every analytic function that is a sixth order or less
polynomial.

For example, choose nx=14 and thus solve 12 equations in 12 unknowns.
Zero based subscripting is chosen for this specific example.

In general the x points are provided xg[0] through xg[nx].

Each equation is based on a unique solution point, i,  with
the right hand side, RHS, being f(xg[i]).

The diagram below represents a row for equation i with
subscript k indicating the column.
The matrix of derivative weights is shown with the
location in the row where the weights are to be entered.

               index k 
          +--+--+--+--+--+--+--+--+--+--+--+--+
row i     | 0| 1| 2| 3| 4| 5| 6| 7| 8| 9|10|11|
          +--+--+--+--+--+--+--+--+--+--+--+--+

               derivative weights, for each derivative
               for each row index i, the value of k, column, is shown
index      0  1  2  3  4  5  6  index j     mid = n/2 = 3
  i   pt +--+--+--+--+--+--+--+
  0    0 |RH| 0| 1| 2| 3| 4| 5|   k = j - 1  (special for RHS j=1,n-1)
         +--+--+--+--+--+--+--+
  1    1 | 0| 1| 2| 3| 4| 5| 6|   k = j + i - pt
         +--+--+--+--+--+--+--+
  2    2 | 0| 1| 2| 3| 4| 5| 6|   k = j + i - pt
         +==+==+==+==+==+==+==+
  3    3 | 0| 1| 2| 3| 4| 5| 6|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  4    3 | 1| 2| 3| 4| 5| 6| 7|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  5    3 | 2| 3| 4| 5| 6| 7| 8|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  6    3 | 3| 4| 5| 6| 7| 8| 9|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  7    3 | 4| 5| 6| 7| 8| 9|10|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  8    3 | 5| 6| 7| 8| 9|10|11|   k = j + i - pt    pt = mid
         +==+==+==+==+==+==+==+
  9    4 | 5| 6| 7| 8| 9|10|11|   k = j + nx - n -2
         +--+--+--+--+--+--+--+
 10    5 | 5| 6| 7| 8| 9|10|11|   k = j + nx - n -2
         +--+--+--+--+--+--+--+
 11    6 | 6| 7| 8| 9|10|11|RH|   k = j + nx - n -2 (special for RHS j=1,n-1)
nx-3  n-1+--+--+--+--+--+--+--+


the a0*U(x) is a special case, because the matrix is the identity matrix.
a0(x) is added to the diagonal element, when i=k

This example yields the following Java language code for fourth
derivative contribution. The matrix must initially be zero and
all terms are summed.
  n = 7;
  nx = 14;
  rhs = nx-1;
  double cuxxxx[][] = new double[n][n]; // initialized by deriv or nuderiv
  double kf[][] = new double[nx-2][nx-1]; // only interior matrix

  for(int i=0; i<nx-2) // walk rows of matrix kf
  {
    row = i
    kf[row][rhs] = f(xg[row]); // right hand side
    if(i==1) // use point, pt=0 for cuxxxx
    {
      pt = 0;
      offs = -1;
      kf[row][rhs] -= a4(xg[row])*cuxxxx[pt][0]*U(xg[0]); // known boundary
      for(int j=1; j<n; j++)
      {
        k = j + offs;
	kf[row][k] += a4(xg[row])*cuxxxx[pt][j]);
      }
    }
    else if(i==nx-3) // use point, pt=n-1
    {
      pt = n-1;
      offs = nx - n - 2;
      kf[row][rhs] -= a4(xg[nx-1])*cuxxxx[pt][n-1]*U(xg[nx-1]); // boundary
      for(int j=1; j<n; j++)
      {
	k = j + offs;
	kf[row][k] += a4(xg[k])*cuxxxx[n-1][j]);
      }
    }
    else // no RHS case, still pick point
    {
      mid = n/2;
      pt  = mid;
      offs = row-pt;
      if(i>nx-n+2)
      {
        pt = i-nx+n+2;
        offs = nx-n-2;
      }
      for(int j=0; j<n; j++) // general case, no RHS
      {
	k = j+offs;
	kf[row][k] += a4(xg[k])*cuxxxx[pt][j];
      }
    }


For the matrix indices starting at 1, changes are:

               index k 
          +--+--+--+--+--+--+--+--+--+--+--+--+
row i     | 1| 2| 3| 4| 5| 6| 7| 8| 9|10|11|12|
          +--+--+--+--+--+--+--+--+--+--+--+--+

               derivative weights, for each derivative
               for each row index i, the value of k, column, is shown
index      0  1  2  3  4  5  6  index j     mid = n/2 = 3
  i   pt +--+--+--+--+--+--+--+
  1    0 |RH| 1| 2| 3| 4| 5| 6|   k = j   (special for RHS j=1,n-1)
         +--+--+--+--+--+--+--+
  2    1 | 1| 2| 3| 4| 5| 6| 7|   k = j + i - pt
         +--+--+--+--+--+--+--+
  3    2 | 1| 2| 3| 4| 6| 6| 7|   k = j + i - pt
         +==+==+==+==+==+==+==+
  4    3 | 1| 2| 3| 4| 5| 6| 7|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  5    3 | 2| 3| 4| 5| 6| 7| 8|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  6    3 | 3| 4| 5| 6| 7| 8| 9|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  7    3 | 4| 5| 6| 7| 8| 9|10|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  8    3 | 5| 6| 7| 8| 9|10|11|   k = j + i - pt    pt = mid
         +--+--+--+--+--+--+--+
  9    3 | 6| 7| 8| 9|10|11|12|   k = j + i - pt    pt = mid
         +==+==+==+==+==+==+==+
 10    4 | 6| 7| 8| 9|10|11|12|   k = j + nx - n -1
         +--+--+--+--+--+--+--+
 11    5 | 6| 7| 8| 9|10|11|12|   k = j + nx - n -1
         +--+--+--+--+--+--+--+
 12    6 | 7| 8| 9|10|11|12|RH|   k = j + nx - n  (special for RHS)
nx-2  n-1+--+--+--+--+--+--+--+


the a0*U(x) is a special case, because the matrix is the identity matrix.
a0(x) is added to the diagonal element, when i=k

This example yields the following Java language code for fourth
derivative contribution. The matrix must initially be zero and
all terms are summed.
  n = 7;
  nx = 14;
  rhs = nx-1;
  double cuxxxx[][] = new double[n][n]; // initialized by deriv or nuderiv
  double kf[][] = new double[nx][nx+1]; // full matrix, only interior used

  for(int i=1; i<nx-1) // walk internal, non boundary rows of matrix kf
  {
    kf[row][rhs] = f(xg[row]); // right hand side
    if(i==1) // use point, pt=0 for cuxxxx
    {
      kf[row][rhs] -= a4(xg[0])*cuxxxx[0][0]*U(xg[0]); // known boundary
      for(int j=0; j<n-1; j++)
      {
        k = j + 1;
	kf[row][k] += a4(xg[rk])*cuxxxx[0][j]);
      }
    }
    else if(i==nx-2) // use point, pt=n-1
    {
      kf[row][rhs] -= a4(xg[nx-1])*cuxxxx[n-1][n-1]*U(xg[nx-1]); // boundary
      for(int j=0; j<n-1; j++)
      {
	k = j + ;
	kf[row][k] += a4(xg[k]*cuxxxx[n-1][j]);
      }
    }
    else // no RHS case, still pick point
    {
      mid = n/2;
      pt  = mid;
      offs = row-pt;
      if(i>nx-n+2)
      {
        pt = i-nx+n+1;
        offs = nx-n-1;
      }
      for(int j=0; j<n; j++) // general case, no RHS
      {
	k = j+offs;
	kf[row][k] += a4(xg[k])*cuxxxx[pt][j];
      }
    }