How to code this iterative scheme for solving 2-group diffusion equation?

I am trying to write a python code to solve the neutron diffusion equation to model neutron flux distribution in a one-dimensional two-group setting. The governing equations of the system are:

$$-D_1\nabla^2\phi_1+(\mathcal{E}_{a_1}+\mathcal{E}_{s12})\phi_1=\frac{1}{k}[\nu_1\mathcal{E}_{f_1}\phi_1+\nu_2\mathcal{E}_{f_2}\phi_2]+\mathcal{E}_{s21}\phi_2$$

$$--D_2\nabla^2\phi_2+(\mathcal{E}_{a_2}+\mathcal{E}_{s21})\phi_2=\mathcal{E}_{s12}\phi_1$$

where $$D$$=Diffusion Co-efficient, $$\phi$$=neutron flux, $$\mathcal{E}_{a}$$=absorption cross-section, $$\mathcal{E}_{s}$$=scattering cross-section, $$\nu$$=neutrons per fission, $$\mathcal{E}_{f}$$=fission cross-section and $$k$$=multiplication factor.

I divided the geometry into 420 mesh elements and discretized the equations using Forward Difference Method (FDM). Putting the whole system into matrix form yields:

$$\Bigl[C_1\Bigr]\Bigl[\phi_1\Bigr]=\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_1\mathcal{E}_{f_1}\Bigr]\Bigl[\phi_1\Bigr]+\frac{1}{k}\Bigl[h\Bigr]\Bigl[\nu_2\mathcal{E}_{f_2}\Bigr]\Bigl[\phi_2\Bigr]+\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s21}\Bigr]\Bigl[\phi_2\Bigr]$$

$$\Bigl[C_2\Bigr]\Bigl[\phi_2\Bigr]=\Bigl[h\Bigr]\Bigl[\mathcal{E}_{s12}\Bigr]\Bigl[\phi_1\Bigr]$$

Where $$\Bigl[C\Bigr]$$=Co-efficient Matrix and $$h$$=mesh element length.

Now, for the $$n$$th iteration, I am supposed to use $$\phi_{1}^{n-1}$$, $$\phi_{2}^{n-1}$$ and $$k^{n-1}$$ to solve the two above equations and calculate $$\phi_{1}^{n}$$, $$\phi_{2}^{n}$$ and $$k^{n}$$ using the iterative scheme:

$$\Bigl[C_1\Bigr]\Bigl[\phi_1\Bigr]^n=\frac{1}{k}\Bigl[\nu_1\mathcal{E}_{f_1}\Bigr]\Bigl[\phi_1\Bigr]^{n-1}+\frac{1}{k}\Bigl[\nu_2\mathcal{E}_{f_2}\Bigr]\Bigl[\phi_2\Bigr]^{n-1}+\Bigl[\mathcal{E}_{s21}\Bigr]\Bigl[\phi_2\Bigr]^{n-1}$$

$$\Bigl[C_2\Bigr]\Bigl[\phi_2\Bigr]^n=\Bigl[\mathcal{E}_{s12}\Bigr]\Bigl[\phi_1\Bigr]^n$$

$$k^n=k^{n-1}\frac{\int dr(\nu_1\mathcal{E}_{f_1}\phi_1^n+\nu_2\mathcal{E}_{f_2}\phi_2^n)}{\int dr(\nu_1\mathcal{E}_{f_1}\phi_1^{n-1}+\nu_2\mathcal{E}_{f_2}\phi_2^{n-1})}$$

I am confused about how to code the last equation. Do I take the differential length $$dr$$ as mesh element length $$h$$ and turn the integration into a summation over the geometry? Or do I actually try some kind of integration? I cannot imagine how an integration can be introduced here given that $$\phi_1$$ and $$\phi_2$$ are not known as functions but as matrices.