Using solution of PDE to find statistical moments

I’m solving a PDE with Dirac delta initial condition

\frac{\partial P}{\partial L}(L,\eta) = \frac{1}{L_{loc}}\bigg(2\eta\frac{\partial P}{\partial \eta} +(\eta^2+1)\frac{\partial^2 P}{\partial \eta^2}\bigg), \quad \eta >1,

with initial condition P(L=0) = \delta(\eta-1). With the help of Jeremy B I have code which solves this

"""
1d Fokker-Planck equation, spatially dependant diffusion process with
Dirac delta initial condition, Dirichlet and Neumann.
"""

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt

L = 2          # final L
num_steps = 10   # number of length steps
dl = L / num_steps # length step size
theta = 1 #Lloc constant term


# Create mesh and define function space
nx = 100
mesh = IntervalMesh(nx, 1., 3)

V = FunctionSpace(mesh, 'P', 1)

x = SpatialCoordinate(mesh)

# Define boundary condition
u_D = Expression("1./(pow(x[0]-y,2) + pow(eps, 2))", eps=1e-6, pi=pi, y=1., degree=2)
# plot(u_D)

def right(x, on_boundary):
    return near(x[0], 1.)
def left(x, on_boundary):
    return near(x[0], 1.)

bc = DirichletBC(V, Constant(0), right)

# Define initial value
u_n = interpolate(u_D, V)
C = assemble(u_n*dx)
u_n.assign(u_n/C)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

a_form = (u*v + (dl/theta)*(x[0]*x[0]-1)*u.dx(0)*v.dx(0))*dx
L_form = u_n*v*dx

# Time-stepping
u = Function(V)
l = 0
for n in range(num_steps):

    # Update current time
    l += dl

    # Compute solution
    solve(a_form == L_form, u, bc)
    
    vtkfile = File('foquesol.pvd')
    vtkfile << u

    # Plot solution
    plot(u)
    #plot(mesh)

    # Update previous solution
    u_n.assign(u)
    print(assemble(u*dx))
   

# Hold plot
plt.xlabel('$\\eta$')
plt.ylabel('$P(\\eta,L)$')
plt.legend(['$P(\\eta,L_0)$','$P(\\eta,L_1)$'], loc='upper right')
plt.show()


I want to use the solution from this code to calculate statistical moments, namely the integral

\int_{1}^{\infty}\left(\frac{2}{1+\eta}\right)^{n} P(L, \eta) d \eta

I know how to do this in python using scipy where I know what the function P is, but not in this case where it is a vector. The resulting plot I want is something like this, where I plot the above integral with varying L on the x axis

Sorry if this is a really dumb question but I’m stuck on the best way to try and achieve this.