Residuals of PDE and numerical vs. symbolic differentiation: |h | 0.1 | 0.01 | 0.001 | 0.0001| |-------------- | ---------------| ---------------- | -----------| ----------------| |residual_u | 0.264407 | 0.000264403 | 4.94253e-07 | 3.97963e-05| |residual_v | 1.76263e-05 | 1.76271e-08 | 3.29687e-11 | 2.55261e-09| |residual_f1 | 36.0141 | 0.00360142 | 3.58131e-07 | 1.55126e-09| |residual_f2 | 0 | 0 | 0 | 0| |residual_u_tld| 37879.2 | 948786 | 1.29167e+06 | 750000| |residual_v_tld| 1.68883 | 1.31748 | 0.319637 | 0.106426| |residual_phi | 34431.6 | 965059 | 1.29569e+06 | 751592| |residual_mu | 1.68883 | 1.31748 | 0.319637 | 0.106429| Obviously, something is wrong with utild and vtild.

(\Delta-\lambda_1^2 I) \int_\Omega G_1(x-y) 1 dy = 1?

Interesting. How much debugging have you done on the volume potential? Could you check whether e.g. ```math (\Delta-\lambda_1^2 I) \int_\Omega G_1(x-y) 1 dy = 1? ```

Barely minimum. I only printed target points and centers etc and they looked fine. Also, I tested on $`f_1(\phi)`$ and $`f_2(\phi)`$ as shown above so the only place that things can go wrong is the volume integral. I will post updates after checking what you suggested.