Hi, I am trying to solve a higher order Eigenproblem, i.e. one which contains fourth order derivatives. Specifically, as a simple example I tried to find the Eigenvalues of the following:

Which should have eigenvalues, \lambda_n = n^2\pi^2 . My attempt to solve this was to introduce a new variable u = v'', and create the weak form from these two equations together. But the answers I am getting are incorrect. Here is a MWE of the code:

```
from fenics import *
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as la
# Define the mesh
Nmesh = 1001
mesh = IntervalMesh(Nmesh, -1.00, 1.00)
# Define the function space
F1 = FiniteElement('P', "interval", 1)
element = MixedElement([F1, F1])
V = FunctionSpace(mesh, element)
# Set Boundary Conditions for v and u
bcs = [DirichletBC(V.sub(0), Constant((0.0)), 'on_boundary'),\
DirichletBC(V.sub(1), Constant((0.0)), 'on_boundary')]
p = TrialFunction(V)
f = TestFunction(V)
# Split functions as necessary
v, u = split(p)
f_1, f_2 = split(f)
# Define the bilinear form
a = (-u.dx(0)*f_1.dx(0) +u*f_2+v.dx(0)*f_2.dx(0))*dx
m = (u*f_1)*dx
# Create Matricies
A, M = PETScMatrix(), PETScMatrix()
b = PETScVector()
# Assemble the matricies
assemble(a, tensor=A)
assemble(m, tensor=M)
# Appply boundary conditions
[bc.apply(M) for bc in bcs]
bcinds = []
for bc in bcs:
bcdict = bc.get_boundary_values()
bcinds.extend(bcdict.keys())
# This just converts PETSc to CSR
A = sp.csr_matrix(A.mat().getValuesCSR()[::-1])
M = sp.csr_matrix(M.mat().getValuesCSR()[::-1])
# Create shift matrix
shift = 1.2345e10*np.ones(len(bcinds))
S = sp.csr_matrix((shift, (bcinds, bcinds)), shape=A.shape)
# Compute Eigenvalues
ev, ew = la.eigs(A+S, 6, M, sigma=0.0,which='SM')
print(ev)
```