Color Mesh Incompatible Dimensions

Question

Error Message

File "/Applications/Spyder.app/Contents/Resources/lib/python3.9/matplotlib/axes/_axes.py", line 5575, in _pcolorargs raise TypeError('Dimensions of C %s are incompatible with'

TypeError: Dimensions of C (101, 101) are incompatible with X (101) and/or Y (101); see help(pcolormesh)

Problem https://www.chegg.com/homework-help/questions-and-answers/problem-1-compute-steady-state-concentration-umol-um3-distribution-channel-shown-figure--p-q96105354

Code

# Import the required modules
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

# Constants
D = 0.5 # (μm^2)/s 
k = 0.1 # 1/s
Lx = 30 # μm
Ly = 10 # μm
Ca = 40 # Initial concentration (μmol/μm^3) 
Cb = 100 # Final concentration (μmol/μm^3) 
Cguess = 70 # Guess concentration (μmol/μm^3)

# Iteration parameters
maxit = 2000
tol = 0.0001 # Relative tolerance
merr = 1e5
lam = 1.4 # Parameter for convergence rate 

# Setup grid
npts = 100
dx = Lx/npts;
xx = np.linspace(0,Lx,num=npts+1,endpoint=True);
nx = npts+1;
dy = Ly/npts;
yy = np.linspace(0,Ly,num=npts+1,endpoint=True);
ny = npts+1;

# Step 1 - Initial Guesses
M = np.ones((ny,nx));    # set matrix size 
M = Cguess*M;    # initial guess for Ca in all nodes

# Step 2 - Apply Boundary Conditions
M[0,:] = Ca
M[:,0] = Cb

# Step 3 and 4 - Apply LDE, walking over nodes
cc = 0 # Counter
a = 1+k*dx*dx/(4*D)
while merr>tol:
    Mold = np.copy(M) # Save current values to old 
    M[-1,-1] = (2*M[-2,-1]+2*M[-1,-2])/(4*a)
    for j in range(1,ny-1):
        for i in range(1,nx-1): 
            M[j,i] = (M[j,i-1]+M[j,i+1]+M[j-1,i]+M[j+1,i])/(4*a); 
    for j in range(1,ny-1):
        M[j,-1] = (2*M[j,-2]+M[j-1,-1]+M[j+1,-1])/(4*a);
    for i in range(1,nx-1): 
        M[-1,i] = (M[-1,i-1]+M[-1,i+1]+2*M[-2,i])/(4*a);
    for i in range(1,68): 
        M[0,i] = (M[0,i-1]+M[0,i+1]+2*M[1,i])/(4*a); 
    M = lam*M+(1-lam)*Mold # Adjust for convergence rate 
    cc = cc+1
    ea = np.abs((M-Mold)/M)
    merr = np.max(ea)
    if cc>maxit: break
avg = np.mean(M[:,0])

# Plot 2D
X,Y = np.meshgrid(xx,yy)
ax = plt.axes(projection="3d")
ax.plot_surface(X,Y,M,edgecolor="none")
ax.view_init(elev=40,azim=-35)
plt.show()

# Plot color mesh
p = plt.pcolormesh(xx, yy, M, cmap="RdBu", shading="flat", vmin=0, vmax=100)
ct = plt.contour(X, Y, M, cmap="gray", levels=10, vmin=0, vmax=100)
c = plt.colorbar(p)
plt.xlabel("Lx (μm)")
plt.ylabel("Ly (μm)")
c.set_label("Concentration (μmol/μm3)")
plt.show()

# Print Results
print("Converged in %d iterations" % cc)
print("Max error is %f" % merr)
print("Mean concentration along central axis = %f μmol/μm^3" % avg)