Mixed Implementation for a 4th order equation

Errors
1

Hello, everyone!
I have a couple PDE who it’s very similar to the Navier-Stokes equations, and I’ve trying to implement a IPCS scheme to my set of equations, however, there is a step in which I have a fourth order operator, and due to that I introduce an auxiliary variable, just like in the following equations:

\frac{1}{Pe}\partial_k^{2}R_i^{n+1} + (\frac{2\lambda}{3} - \frac{\alpha}{\Delta t})R_{i}^{n+1} = \partial_j(\frac{1}{Pe}\partial_k^{2}Q_{ij}^{\star} - \frac{1}{\Delta t}Q_{ij}^{\star}) + \frac{2\lambda}{3} R_i^{n},

And we have the following relationship between R_i and u_i, which reads

R_i^{n+1} = \frac23 \partial_j^2u_i^{n+1}

Essentially, the right-hand side of the following equation is known, because I start the simulation with a given flow, and Q_{ij}^{\star} is a quantity calculated before.
I have been trying to implement something like in the Mixed Poisson Demo, for the previous set of equations, considering the contribution from Q^{\star} = 0 , just for simplicity, but I get the following error ValueError: not enough values to unpack (expected 1, got 0)
Here’s my code:

## Define parameters for time stepping and constants
dt = 0.1
num_steps = 1
T = dt *num_steps            ## final time
lamb = 0.1         ## lambda parameter
Pe   = 10.0         ## Peclet number

## Create the mesh
x_channel = 0.5
y_channel = 0.2
nx = ny = 20
mesh = RectangleMesh(Point(0, 0), Point(x_channel, y_channel), 2*nx, ny, 'crossed')

## Define vector elements and function spaces
UE = VectorElement("Lagrange", mesh.ufl_cell(), 2)

# Define function spaces
ME = MixedElement([UE, UE])
U  = FunctionSpace(mesh, ME)

Rspace = U.sub(0)
Uspace = U.sub(1)

## Define boundaries
inflow   = 'near(x[0], 0)'                     
outflow  = 'near(x[0], 0.5)'
walls    = 'near(x[1], 0) || near(x[1], 0.2)'

## Define inflow profile for the velocity
inflow_profile = ('4.0*1.5*x[1]*(0.1 - x[1]) / pow(0.1, 2)', 
                  '0')

## Define boundary conditions for the velocity
## For mixed systems, we should specify which subsystem we set
## the boundary condition. This is done via the sub(0) and sub(1)
bcu_inflow = DirichletBC(Uspace, 
                         Expression(inflow_profile, 
                                    degree=2), 
                         inflow)
bcu_walls = DirichletBC(Uspace, 
                        Constant((0, 0)), 
                        walls)
bcu = [bcu_inflow, bcu_walls]         ## Colect all boundary conditions for the velocity

## Define trial and test functions
(R, u)   = TrialFunctions(U)
(tau, v) = TestFunctions(U)

## Define expressions used in variational forms
k1 = Constant(dt)
k2 = Constant(Pe)

F = (1 / k2) * inner(nabla_grad(R), nabla_grad(tau))*dx \
   + (2*lamb / 3. - (1 / k1)) *  inner(R, tau)*dx \
   + (2 / 3) * inner(nabla_grad(u), nabla_grad(v))*dx - inner(R, v)*dx

a2 = lhs(F)
L2 = rhs(F)

# Assemble matrices
A2 = assemble(a2)

# Apply boundary conditions to matrices
[bc.apply(A2) for bc in bcu]

# Time-stepping
w = Function(U)

t = 0
for n in range(num_steps):
    
    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b2 = assemble(L2)
    [bc.apply(A2) for bc in bcu]
    solve(A2, w, b2)
    (R, u) = w.split()

Hope you could help me.
Thanks