A third order equation of the form:

and

I haven't tried this, but this should be possible with a coupled equation in a fully implicit manner. Here is an attempt. The sign on the diffusion term seems to make a large difference in the stability.

from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, Viewer
from fipy import UpwindConvectionTerm, ImplicitSourceTerm

m = Grid1D(nx=100, Lx=1.)

phi = CellVariable(mesh=m, hasOld=True, value=0.5, name=r'$\phi$')
psi = CellVariable(mesh=m, hasOld=True, name=r'$\psi$')

phi.constrain(0, m.facesLeft)
phi.constrain(1, m.facesRight)

psi.constrain(0, m.facesRight)
psi.constrain(0, m.facesLeft)

eqnphi = TransientTerm(1, var=phi) == UpwindConvectionTerm([[1]], var=psi)
eqnpsi = ImplicitSourceTerm(1, var=psi) == -DiffusionTerm(1, var=phi)

eqn = eqnphi & eqnpsi

vi = Viewer((phi, psi))

for t in range(1000): 
    phi.updateOld()
    eqn.solve(dt=0.001)
    print phi
    vi.plot()


raw_input('stopped')