# Stiffness Matrix Formation for PDE with Neumann Boundary

Given the problem $$-\nabla u + u = f$$ $$n\cdot\nabla u = g \quad\text{on} \quad \Gamma$$ I can show the discretization given through the Galerkin formulation is $$Au=b$$ where $$A = \int_\Omega \nabla\varphi_i\cdot\nabla\varphi_jdx + \int_\Omega \varphi_i\varphi_jdx$$ $$b= \int_\Omega f \varphi_idx + \int_\Gamma g\varphi_ids$$ and the basis functions are linear with form $$\varphi_i = c_1+c_2x+c_3y$$, satisfying the condition that $$\varphi_i(x_j) = \delta_{ij}$$ For $$A$$, I can find the first term's elemental stiffness matrix through $$A = \int_\Omega (c_{2,i}c_{2,j} + c_{3,i}c_{3,j})dx = Area(c_{2,i}c_{2,j} + c_{3,i}c_{3,j})$$

I'm trying to calculate the elemental stiffness matrix for $$\int_\Omega\varphi_i\varphi_jdx$$ but getting stuck as I'm not sure what do with $$c_3$$ terms

Since $$\varphi_i \varphi_j$$ is a polynomial of degree 2, you can just use a Gauss quadrature formula (if the formula you choose has degree $$\ge 2$$, it will yield the exact value for the integral). You can either compute the quadrature nodes form each element, or compute all the integrals in a reference element, using a change of variables.