Summary: We study a system of parabolic equations with nonlocal nonlinearity of the following type \[ \begin{cases} u_t - a_1(l_1(u), l_2(v)) \Delta u + \lambda_1 | u |^{p - 2} u = f_1(x, t) & \text{in } \Omega \times] 0, T] \\ v_t - a_2(l_1(u), l_2(v)) \Delta v + \lambda_2 | v |^{p - 2} v = f_2(x, t) & \text{in } \Omega \times] 0, T] \\ u(x, t) = v(x, t) = 0 & \text{on } \partial \Omega \times] 0, T] \\ u(x, 0) = u_0(x), \, v(x, 0) = v_0(x) & \text{in } \Omega, \end{cases} \] where \(a_1\) and \(a_2\) are Lipschitz-continuous positive functions, \(l_1\) and \(l_2\) are continuous linear forms, \(\lambda_1, \lambda_2 \geq 0\) and \(p \geq 2\).
We prove the convergence of a linearized Euler-Galerkin finite element method and obtain the order of convergence in the \(L_2\) norm. Finally we implement and simulate the presented method in Matlab’s environment.