import numpy as np

nx=3
nz=4
flux2 = np.empty((nx+1, nz+1))
flux2[:] = np.random.normal(size=(nx+1, nz+1))
flux3 = flux2.copy()
bz = np.random.normal(size=(nx+1, nz+1))
bx = np.random.normal(size=(nx+1, nz+1))

dz = 0.112344
dx = 0.015543
for i in np.arange(1,nx+1):
    for j in np.arange(1,nz+1):
        flux2[i,j] = 0.5*(flux2[i-1,j] + bz[i-1,j]*dx) + 0.5*(flux2[i,j-1] - bx[i,j-1]*dz)

a = 0.5
aexp = np.arange(nz).reshape(nz, 1) - np.arange(nz).reshape(1, nz)
abcoeff = 0.5**aexp
abcoeff[aexp<0] = 0
for i in np.arange(1,nx+1):
    b = 0.5*flux3[i-1, 1:] + 0.5*bz[i-1, 1:]*dx - 0.5*bx[i,:-1]*dz
    bvals = (abcoeff * b.reshape(1, nz)).sum(axis=1)
    n = np.arange(1, nz+1)
    x0 = flux3[i, 0]
    #ab = b * np.cumsum(a**(n-1))
    flux3[i, 1:] = a**n * x0 + bvals

print np.abs(flux2-flux3) < 1e-15