Sim.plot() doesn't seem to work for 1D variables of the BaseModel Class

Hi, I’ve been trying to experiment with building models from scratch using the PyBaMM.BaseModel class. However, I am unable to plot 1D variables this way.

My model is shared below here:

import pybamm 

class model(pybamm.BaseModel):
    
    def __init__(self):
        super().__init__()
        pybamm.citations.register("Marquis2019")
        self.rhs = {}
        x = pybamm.Variable("x", domain="negative electrode")
        y = pybamm.Variable("y")
        z = pybamm.Variable("z")
        x_surf = pybamm.surf(x)
        z_eq = pybamm.Scalar(3) * x_surf + pybamm.Scalar(4) * y +z
        self.algebraic[z] = z_eq
        a = pybamm.Parameter("MyParam")
        dxdt = pybamm.Scalar(1)*x + y
        dydt = y
        self.boundary_conditions = {x: {
            "left": (pybamm.Scalar(0), "Dirichlet"),
            "right": (pybamm.Scalar(0), "Neumann"),
        }}
        # define geometry
        xs = pybamm.SpatialVariable("X_s", domain=["negative electrode"], 
            coord_sys="cartesian"
        )
        # mesh and discretise
        submesh_types = {"negative electrode": pybamm.Uniform1DSubMesh}
        var_pts = {xs: 20}
        
        geometry = {"negative electrode": {xs: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
        spatial_methods = {"negative electrode": pybamm.FiniteVolume()}
        disc = pybamm.Discretisation(mesh, spatial_methods)


        self.rhs = {x: dxdt, y: dydt}
        self.initial_conditions = {x: pybamm.Scalar(1), y: pybamm.Scalar(2), z: pybamm.Scalar(0)}
        self.variables = {"x": x, "y": y, "z": z}
        self = disc.process_model(self)

And this is the code I am using to run the model:

from MyModel2 import model
import pybamm
from parameter_values import get_parameter_values
import matplotlib
matplotlib.use("TkAgg")  # Use TkAgg backend for interactive plotting


model1 = model()
# parameter_values = pybamm.ParameterValues("Chen2020")
# parameter_values.update({"R_c": 2, "C_c": 60}, check_already_exists=False)
#parameter_values = pybamm.ParameterValues(get_parameter_values())
parameter_values = pybamm.ParameterValues({"MyParam": 1})
sim = pybamm.Simulation(model1, parameter_values=parameter_values, solver = pybamm.IDAKLUSolver(), 
                        spatial_methods={"negative electrode" : pybamm.Uniform1DSubMesh})
solution=sim.solve([0, 5])

sim.plot(["x"]) 
sim.plot(["y"])

I am getting an error while plotting “x”, which is a 1D variable. Does anyone have advice on how I can solve this?