Why is facecolor argument in plot_surface() of matplotlib not working in python?

I want my python code to generate a plot with alternating colors. For the sake of the problem, the alternating colors are irrelevant so let's say they're red and blue. I'm creating a numpy array of the alternating colors with the same shape as the inputs x, y and z (251, 251)

colortuple = ('red', 'blue')
colors = np.array(
    [np.array([colortuple[int((xn + yn) % len(colortuple))] for xn in range(xrange)]) for yn in range(trange)])

The array prints the following to the console:

[['red' 'blue' 'red' ... 'red' 'blue' 'red']
 ['blue' 'red' 'blue' ... 'blue' 'red' 'blue']
 ['red' 'blue' 'red' ... 'red' 'blue' 'red']
 ...
 ['red' 'blue' 'red' ... 'red' 'blue' 'red']
 ['blue' 'red' 'blue' ... 'blue' 'red' 'blue']
 ['red' 'blue' 'red' ... 'red' 'blue' 'red']]

which is what I want.

However, when I pass it to plot_surface() as the facecolors argument

axes.plot_surface(X=x, Y=y, Z=z, facecolors=colors)

the plot displayed only has the color of the first argument in colortuple, in this case red.

All 4 variables {x, y, z, colors} have the same (251, 251) shape and are all numpy arrays. I motivated my code from the following website: https://matplotlib.org/2.0.2/mpl_toolkits/mplot3d/tutorial.html

The problem is similar to facecolors in plot_surface matplotlib however I don't get any errors in the console, just a faulty execution.

What am I doing wrong?

Edit:

When I first posted this issue I was too inexperienced to know of the significance of reproducibility. In case anybody else would like to reproduce this issue, here is a slightly modified version of the original code:

from matplotlib.ticker import LinearLocator
from matplotlib import pyplot as plt
from numpy import pi, sin, array, meshgrid, zeros, linspace

# Constants
xmax = 250
xrange = range(xmax + 1)
c = 20.0
dx = 1 / xmax
dt = 0.05 / xmax
k = 4 * pi
multiplier = (c * c * dt * dt / (dx * dx))

# Position and Time Steps.
xsteps = linspace(0, 1, xmax + 1)
tsteps = linspace(0, xmax * dt, xmax + 1)

# Function Psi
# For explanation look up the manual calculations.
psi = zeros(shape=(len(xrange), len(xrange)))
psi[0] = array([
    0 if (xn == 0 or xn == xmax) else
    5 * sin(k * xsteps[xn]) for xn in xrange
])
psi[1] = array([
    0 if (xn == 0 or xn == xmax) else
    5 * sin(k * xsteps[xn]) + multiplier
    * (2.5 * sin(k * xsteps[xn - 1]) + 2.5 * sin(k * xsteps[xn + 1]) - 5 * sin(k * xsteps[xn]))
    for xn in xrange
])
for tn in range(1, xmax):
    psi[tn + 1] = array([
        0 if (xn == 0 or xn == xmax)
        else 2 * psi[tn, xn] - psi[tn - 1, xn] + multiplier * (psi[tn, xn - 1] + psi[tn, xn + 1] - 2 * psi[tn, xn])
        for xn in xrange
    ])

# Setting the parameters for plotting via matplotlib
axes = plt.figure().add_subplot(projection="3d")
(x, y) = meshgrid(xsteps, tsteps)

# Create an empty array of strings with the same shape as the meshgrid,
# and populate it with two colors in a checkerboard pattern.
colortuple = ("red", "blue")
colors = array([array([colortuple[int((xn + yn) % len(colortuple))] for xn in xrange]) for yn in xrange])

# Plotting function.
axes.plot_surface(x, y, psi, facecolors=colors)

# Customize the z axis.
axes.zaxis.set_major_locator(LinearLocator(6))
axes.set_xlabel(r"$x\,[\mathrm{m}]$")
axes.set_ylabel(r"$t\,[\mathrm{s}]$")
axes.set_zlabel(r"$\psi\,(x,t)$")
plt.tight_layout()
plt.show()

The problem is that if you don't set the rcount and ccount parameters the renderer automatically downgrades the number of patches to a default maximum of 50 rows/columns, and the downsampling process can pick colors somewhat unpredictably. If you do set them to your actual color array dimensions, the render is slow, beautiful in high res and correct.

From https://matplotlib.org/stable/api/_as_gen/mpl_toolkits.mplot3d.axes3d.Axes3D.plot_surface.html :

The rcount and ccount kwargs, which both default to 50, determine the maximum number of samples used in each direction. If the input data is larger, it will be downsampled (by slicing) to these numbers of points.