Formulation weak form of PDE for solving in FeniCS

I'm trying to derive the weak form for the following Fokker-Planck equation to use in a finite element package. I have

$$\frac{\partial P}{\partial L}(L,\eta) = \frac{1}{\theta}\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)$$. The Dirichlet boundary conditions can be set later. I require the problem to be in the weak form (page 11) for it to be solved by the FeniCS finite element package. I choose to discretize the time derivative via

$$\bigg(\frac{\partial P}{\partial L}\bigg)^{n+1} \approx \frac{P^{n+1}-p^n}{\Delta L}$$

so

$$P^{n+1}-\frac{\Delta L}{\theta}\bigg(2\eta\frac{\partial P^{n+1}}{\partial\eta}+(\eta^2+1)\frac{\partial^2 P^{n+1}}{\partial \eta^2}\bigg) - p^n = 0.$$

The weak form of the problem I want to write like

$$a(P,v) = L_{n+1}(v)$$

Where $$v$$ is a test function. I think then I can write

\begin{align} a(P,v) &= \int_{\Omega}P^{n+1}v-\frac{\Delta L}{\theta}\bigg(2\eta\frac{\partial P^{n+1}}{\partial\eta}+(\eta^2+1)\frac{\partial^2 P^{n+1}}{\partial \eta^2}\bigg)\,dx, \\ L_{n+1}(v) &= \int_{\Omega}p^{n}v\,dx. \end{align} I'm not sure whether this is correct, and if I should do something with the second derivative in $$a(P,v).$$

Actually, the boundary conditions may change the expression of the weak form. Here I am going to assume that you have homogeneous Dirichlet boundary conditions at $$\eta=1$$ as well as at $$\eta_\max>1$$ which is where your are going to truncate your (at the moment unbounded, so untractable numerically) domain in $$\eta$$.
Observe that your equation rewrites $$\frac{\partial P}{\partial L}(L,\eta) = \frac{1}{\theta}\frac{\partial}{\partial \eta}\bigg( (\eta^2+1)\frac{\partial P}{\partial \eta}\bigg)$$
So after your backward Euler discretization in time, it yields $$P^{n+1} - \frac{\Delta L}{\theta}\frac{\partial}{\partial \eta}\bigg( (\eta^2+1)\frac{\partial P^{n+1}}{\partial \eta}\bigg) = P^n$$ in strong form.
To obtain the weak form you can formally integrate the equation against a test function $$v$$ that vanishes at the endpoints of your interval. It yields after one integration by parts in the second term of the left-hand side
$$\int_1^{\eta_\max}P^{n+1}(\eta)v(\eta)d\eta + \int_1^{\eta_\max}\frac{\Delta L}{\theta}(\eta^2+1)\frac{\partial P^{n+1}}{\partial \eta}(\eta)\frac{\partial v}{\partial \eta}(\eta)d\eta = \int_1^{\eta_\max}P^n(\eta)v(\eta)d\eta$$
From there you have your bilinear form $$a(P,v)=\int_1^{\eta_\max}P(\eta)v(\eta)d\eta + \int_1^{\eta_\max}\frac{\Delta L}{\theta}(\eta^2+1)\frac{\partial P}{\partial \eta}(\eta)\frac{\partial v}{\partial \eta}(\eta)d\eta$$ and your linear form $$L_n(v)=\int_1^{\eta_\max}P^n(\eta)v(\eta)d\eta$$