How to solve this time-dependent PDE in Python?

Is there a way to numerically solve the following PDE in Python?

The second term on the RHS has a derivative with respect to time as well as space.

I tried using Py-PDE package in Python, it solves only the form dy(x,t)/dt = f(y(x,t)) so I tried to use a root finding algorithm similar to scipy fzero to get the solution to dy(x,t)/dt - f(y(x,t),dy(x,t)/dt) = 0 (solving for dy(x,t)/dt).

class nonlinearPDE(pde.PDEBase):
    def __init__(self, bc={"derivative":0}):
        self.bc = bc #boundary conditions for operators
    
    def _make_pde_rhs_numba(self, state):
        """numba-compiled implementation of the PDE"""
        laplace = state.grid.make_operator("laplace", bc=self.bc)
        def findroot(f, df, x0, nmax):
            """Newton–Raphson method"""
            for i in range(nmax):
                x0 = x0 - f(x0)/df(x0)
            return x0
    
        @jit
        def pde_rhs(y, t):
            func = lambda dydt : dydt - a*laplace(y) - b*laplace(dydt)
            dydt = findroot(func, lambda x : 1, 0, 1)
            return dydt
    return pde_rhs

However, when the program tries to solve the PDE it throws an error:

  File "...\anaconda3\lib\site-packages\pde\solvers\controller.py", line 191, in run
    t = stepper(state, t, t_break)

  File "...\anaconda3\lib\site-packages\pde\solvers\scipy.py", line 82, in stepper
    sol = integrate.solve_ivp(

  File "...\anaconda3\lib\site-packages\scipy\integrate\_ivp\ivp.py", line 542, in solve_ivp
    solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)

  File "...\anaconda3\lib\site-packages\scipy\integrate\_ivp\rk.py", line 94, in __init__
    self.f = self.fun(self.t, self.y)

  File "...\anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 138, in fun
    return self.fun_single(t, y)

  File "...\anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 20, in fun_wrapped
    return np.asarray(fun(t, y), dtype=dtype)

  File "...\anaconda3\lib\site-packages\pde\solvers\scipy.py", line 74, in rhs_helper
    return rhs(state_flat.reshape(shape), t).flat  # type: ignore

  File "...\anaconda3\lib\site-packages\numba\core\dispatcher.py", line 420, in _compile_for_args
    error_rewrite(e, 'typing')

  File "...\anaconda3\lib\site-packages\numba\core\dispatcher.py", line 361, in error_rewrite
    raise e.with_traceback(None)

TypingError: Cannot capture the non-constant value associated with variable 'y' in a function that will escape.

Since no one has posted an answer yet, I managed to get a minimal working example by using scipy odeint with a method of lines to solve the PDE, that is, by discretizing the Laplace operator, and then wrapping the differential equation inside fsolve to get dydt:

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import fsolve

a=1
b=1
t = np.linspace(0,1,100)
y0 = np.asarray([1,0.5,0]) #starting values
N=len(y0)

def laplace(y):
    """ Laplace operator in cartesian coordinate system """
    d2y_0 = [2*(y[1] - y[0])] #at i=0
    d2y_N = [2*(y[-2] - y[-1])] #at i=N
    d2y_i = [y[i+2] - 2*y[i+1] + y[i] for i in range(N-2)] #elsewhere
    return np.asarray(d2y_0 + d2y_i + d2y_N)

def rhs(y, dydt, t):
    """ RHS of PDE including dydt term """
    return a*laplace(y) + b*laplace(dydt)

def pde(y, t):
    """ Function that solves for dydt using root finding """
    return fsolve(lambda dydt : dydt - rhs(y, dydt, t),np.zeros(N))

#calculate
sol = odeint(pde, y0, t)

#plot odeint with fsolve
import matplotlib.pyplot as plt
plt.figure()
plt.plot(t,sol)
plt.grid(ls='--')
plt.xlabel('$t$')
plt.ylabel('$y_i$')
plt.title('odeint with fsolve')
plt.legend(["$i=$"+str(i) for i in range(N)])

#plot theoretical
u10=y0[0]
u20=y0[1]
u30=y0[2]

E1=np.exp(-2*a*t/(1+2*b)) #1st exponential term
E2=np.exp(-4*a*t/(1+4*b)) #2nd exponential term
u1 = 1/4*((1+2*E1+E2)*u10 - 2*(E2-1)*u20 + (1-2*E1+E2)*u30)
u2 = 1/4*((1-E2)*u10      + 2*(E2+1)*u20 + (1-E2)*u30)
u3 = 1/4*((1-2*E1+E2)*u10 - 2*(E2-1)*u20 + (1+2*E1+E2)*u30)

plt.figure()
plt.plot(t,u1,t,u2,t,u3)
plt.grid(ls='--')
plt.xlabel('$t$')
plt.ylabel('$y_i$')
plt.title('Theoretical')
plt.legend(["$i=$"+str(i) for i in range(N)])

Note that the same Laplace discretization method allows us to solve a system of ODEs to get the exact analytical solution with which we verify our numerical calculation (for a N=3 size system):

>>> np.allclose(sol[:,0],u1)
True
>>> np.allclose(sol[:,1],u2)
True
>>> np.allclose(sol[:,2],u3)
True

Seems it works as intended.