// Integration of (electric field * normal) along line in the distance dint from x-axis 
// solution found with FEM is named u

real qq = 0.0;

real dint = 0.20;
int  nint = 4000;
real h = (a - 2 * dint) / nint;
for(int i = 1; i <= nint; i++) {
  qq = qq + dy(u)(dint + (i - 0.5) * h, dint);
}

qq = h * qq;

cout <<" Charge: "<< qq <<endl;