Tetra-series

[Gottfried]

\( F = \mathbf{B}^{-1}G \)
and when the matrix size is even I get the first series, and when the matrix size is odd, I get the second series.

Hmm, I didn't catch the actual computation, but that's possibly not yet required.
I guess, you need the technique of divergent-summation ; it may be, that the implicite series, which occur by multiplication of \( \mathbf{B}^{-1}G \), have alternating signs, are not converging well, or even are divergent. So you could include the process of Euler-summation.

I found a very nice method to get at least an overview, whether such matrix-product suffers from non-convergence (but which is reparable by Cesaro/Euler-summation).
I've defined a diagonal-matrix dE(o) of coefficients for Euler-summation of order o, which can simply be included in the matrix-product.
Write \( \mathbf{B}^{-1} * dE(o) * G \) where o=1.0 o=1.5 or o=2. With too small o the implicite sums in matrix-multiplication begin to oscillate from some terms (order is "too weak") , for too high o the oscillation is so heavily suppressed, that with dim-number of terms the series is not yet converging. In a matrix-product using dimension dim the number of dim^2 of such sums occur. While likely not all that sums can be handled correctly by that same Euler-vector, for some of them you will see a well approximated result and a general smoothing, making the result matrix-size independent, especially averages between size dim and size dim+1.

In my implementation order o=1 means no change; simply dE(1) = I ; dE(1.7)..dE(2.0) can sum 1-1+1-1... and similar, dE(2.5)..dE(3.0) can sum divergence of type 1-2+4-8+...-... and so on
(details for toy-implementation in Pari/GP see below, I can send you some example scripts if this is more convenient)

Gottfried

Code:

//  Pari/GP
\\ a vector of length dim returning coefficients for Euler-summation of
\\ order "order" (E(1) gives the unit-vector:direct summation
{E(order, dim=n) = local(Eu);
Eu=vector(dim);Eu[1]=order^(dim-1);
for(k=2,dim,Eu[k]=Eu[k-1]-(order-1)^(dim-k+1)*binomial(dim-1,k-2));
Eu=Eu/order^(dim-1);
return(Eu);}

\\ returns this as diagonal-matrix
dE(order,dim=n) = return( matdiagonal(E(order,dim)) )