Runtime error "Unable to define linear variational problem"

Question

I keep getting the following runtime error:

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.
*** Reason:  Expecting the left-hand side to be a bilinear form (not rank 1).
*** Where:   This error was encountered inside LinearVariationalProblem.cpp.
*** Process: 0
*** 
*** DOLFIN version: 1.6.0
*** Git changeset:  dcc25c029d66c32bb714a01bfc832b34c7120e72
*** -------------------------------------------------------------------------

I am using the following code:

    import random
from dolfin import *


# Class representing the intial conditions
class InitialConditions(Expression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(mpi_comm_world()))
    def eval(self, values, x):
        values[0] = 0.63 + 0.02*(0.5 - random.random())
        values[1] = 0.0
    def value_shape(self):
        return (2,)

# Mesh 
mesh = UnitSquareMesh(96,96)

# Mixed function space
U = FunctionSpace(mesh,"Lagrange",1)
W = FunctionSpace(mesh,"Lagrange",1)
V = MixedFunctionSpace([U,W])

# Time variables
dt = 1.0e-1 # time-step
t = float(dt)
T = 2.0 # simulation time

# Current solution
w1 = Function(V)
(c1,mu1) = split(w1)

# Previous solution
w0 = Function(V)
(c0,mu0) = split(w0)

# Variational problem at each time step
(phi,psi) = TestFunctions(V)
w = Function(V)
(c,mu) = split(w)

# Model parameters
M = Constant(1.0)
theta = 0.5
lam = 1.0e-2

# f and df/dc
c = variable(c)
f = 100*c**2*(1-c)**2
df = diff(f, c)

# Create intial conditions and interpolate
w_init = InitialConditions(degree=1)
w0.interpolate(w_init)

# Bilinear form
b1 = c*phi*dx + theta*dt*M*inner(grad(mu),grad(phi))*dx 
b2 = mu*psi*dx - lam*inner(grad(c),grad(psi))*dx
b = b1 + b2
# Linear form
L = c0*phi*dx + (theta - 1.0)*dt*M*inner(grad(mu0),grad(phi))*dx + df*psi*dx

while (t <= T):
    solve(b == L,w1)
    w0.assign(w1)
    t += float(dt)