Full Matrix

It is easy enough to build the marix for the FDM problem.  You just build an Identity matrix of the appropriate size and then replace the rows corresponding to internal points with the appropriate "111-6111" row.  Here it is in action.

Size  =  7;
Δx  =  
        (1. /(Size - 1));
A =  IdentityMatrix[Size^3];
f = Table[ 0,  {(Size^3)}];
Do[
n  =  
        1 + i  +  Size j  +  
          Size^2 k; 
f[[n]] =  Δ(x^2) (i Δx  +  
              j Δx  +  (((k Δx))^2)); 
A[[
          n, {n - (Size^2), n - Size, n - 1, n, n + 1, n + Size, 
            n + (Size^2)}]] = {1. , 1. , 1. , -6. , 
          1. , 1. , 1. }, 
{i, 1, 
        Size - 2}, 
{j, 1, 
        Size - 2}, 
{k, 1, 
        Size - 2}
]

The resulting matrix is very sparse

ListDensityPlot[-Sign[Abs[A]], 
Mesh -> 
          False, 
PlotLabel -> "FDM Matrix"];

Solving the system is simple, but reassembling the solution is slightly more complicated.

VecSol  = LinearSolve[ A,  f];
BlockSol  =  
        Partition[
Partition[
VecSol, 
Size], 
Size];
Dimensions[BlockSol]

Of course, as a prelimnary check if the soln is meaningful we should plot it.

\(Table[
ListPlot3D[\(-\ 
          Sol[[
            k]]\), 
PlotRange -> {0, \  .1}, 
PlotLabel -> 
StringForm["z=``", \ \((k - 
                  1)\)\ Δx]
], 
{k, \ 1, \ Size}]\)

Sparse Matrix

In some ways the sparse version is easier.  In other ways, the book keeping is slightly more complex because we need to work out which points are boundary points.

BoundaryTest[{i_, j_, k_}, Size_]  :=  
      Module[
{}, 
Intersection[{i, j, 
              k}, {0, Size - 
                1}] ≠ {}
]

Once you have done this it is straightforward to build the sparse matrix.

Size  =  15;
Δx  =  
        (1. /(Size - 1));
f = Table[ 0,  {(Size^3)}];
SparseA  = 
        Flatten[
Table[
n  =  
              1 + i  +  Size j  +  
                Size^2 k; 
f[[n]] =  Δ(x^2) (i Δx  +  
                    j Δx  +  (((k Δx))^2)); 
If[
BoundaryTest[{i, j, k}, 
                Size], 
{n, n} -> 
                1, 
{
{n, 
                    n - (Size^2)} -> 
                  1.0, 
{n, n - Size} -> 
                  1.0, 
{n, n - 1} -> 
                  1.0, 
{n, 
                    n} -> -6.0, 
{n, 
                    n + 1} -> 
                  1.0, 
{n, n + Size} -> 
                  1.0, 
{n, n + (Size^2)} -> 
                  1.0
}
], 
{i, 0, Size - 1}, 
{j, 0, 
              Size - 1}, 
{k, 0, 
              Size - 1}
]
];

Note:  this uses a lot less memory than the full matrix for large problems. In addition, it is also faster to solve.

in  =  Date[];
VecSol = 
        Developer`SparseLinearSolve[SparseA, f];
Date[] - in

Of course, we still need to reassemble the solution in a meaningful manner and plot the result to check we did not mess up somewhere.

Length[VecSol]