Solving integral within PDEs and ODE system

Hello all,

I have a 2D dynamic simulation with 4 PDEs. In addition to the PDEs, I must also incorporate an ODE, as follows:

\frac{\mathrm{d} u}{\mathrm{d} t}= \alpha (u_{o}-u)+\beta \int_{0}^{L}[u-u_{w}] \mathrm{d} z

My first approach consisted of compiling the whole Variational Formulation (VF) inside the For loop so I could solve the integral in each iteration (using scipy). However, this coarse solution makes the computation to last at least 20 more times than the case in which I compile the VF outside the loop. I do not know how to set up the variational problem so I could solve the complete system (4 PDEs + 1 ODE), including the integral calculation.

I really appreciate your help on this topic. Any suggestion is welcome.

Best regards,