def METHOD_NAME(plotIt=True, n=60):
    np.random.seed(5)
    # Here we are going to rearrange the equations:
    # (phi_ - phi)/dt = A*(d2fdphi2*(phi_ - phi) + dfdphi - L*phi_)
    # (phi_ - phi)/dt = A*(d2fdphi2*phi_ - d2fdphi2*phi + dfdphi - L*phi_)
    # (phi_ - phi)/dt = A*d2fdphi2*phi_ + A*( - d2fdphi2*phi + dfdphi - L*phi_)
    # phi_ - phi = dt*A*d2fdphi2*phi_ + dt*A*(- d2fdphi2*phi + dfdphi - L*phi_)
    # phi_ - dt*A*d2fdphi2 * phi_ =  dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi
    # (I - dt*A*d2fdphi2) * phi_ =  dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi
    # (I - dt*A*d2fdphi2) * phi_ =  dt*A*dfdphi - dt*A*d2fdphi2*phi - dt*A*L*phi_ + phi
    # (dt*A*d2fdphi2 - I) * phi_ =  dt*A*d2fdphi2*phi + dt*A*L*phi_ - phi - dt*A*dfdphi
    # (dt*A*d2fdphi2 - I - dt*A*L) * phi_ =  (dt*A*d2fdphi2 - I)*phi - dt*A*dfdphi
    h = [(0.25, n)]
    M = discretize.TensorMesh([h, h])
    # Constants
    D = a = epsilon = 1.0
    I = discretize.utils.speye(M.nC)
    # Operators
    A = D * M.face_divergence * M.cell_gradient
    L = epsilon**2 * M.face_divergence * M.cell_gradient
    duration = 75
    elapsed = 0.0
    dexp = -5
    phi = np.random.normal(loc=0.5, scale=0.01, size=M.nC)
    ii, jj = 0, 0
    PHIS = []
    capture = np.logspace(-1, np.log10(duration), 8)
    while elapsed < duration:
        dt = min(100, np.exp(dexp))
        elapsed += dt
        dexp += 0.05
        dfdphi = a**2 * 2 * phi * (1 - phi) * (1 - 2 * phi)
        d2fdphi2 = discretize.utils.sdiag(a**2 * 2 * (1 - 6 * phi * (1 - phi)))
        MAT = dt * A * d2fdphi2 - I - dt * A * L
        rhs = (dt * A * d2fdphi2 - I) * phi - dt * A * dfdphi
        phi = Solver(MAT) * rhs
        if elapsed > capture[jj]:
            PHIS += [(elapsed, phi.copy())]
            jj += 1
        if ii % 10 == 0:
            print(ii, elapsed)
        ii += 1
    if plotIt:
        fig, axes = plt.subplots(2, 4, figsize=(14, 6))
        axes = np.array(axes).flatten().tolist()
        for ii, ax in zip(np.linspace(0, len(PHIS) - 1, len(axes)), axes):
            ii = int(ii)
            M.plot_image(PHIS[ii][1], ax=ax)
            ax.axis("off")
            ax.set_title("Elapsed Time: {0:4.1f}".format(PHIS[ii][0]))