diff --git a/examples/douglas_2024/README.md b/examples/douglas_2024/README.md
index b66c9fa..dc9076d 100644
--- a/examples/douglas_2024/README.md
+++ b/examples/douglas_2024/README.md
@@ -58,9 +58,9 @@ Note: this example does not make use of adaptive meshing.
 ```sh
 FreeFem++-mpi -v 0 examples/douglas_2024/grabowski.md -mo $workdir/G
 ```
-## Run the standalone serial FreeFEM code provided in the Supplementary Materials
+## Run the code provided in the Supplementary Materials
 ```sh
-FreeFem++ -v 0 examples/douglas_2024/example1_suppmat.md
+FreeFem++-mpi -v 0 examples/douglas_2024/example1_suppmat.md
 ```
 
 ## Perform parallel computations for Grabowski--Berger vortex flow using `ff-bifbox`
diff --git a/examples/douglas_2024/basecomputeaug.md b/examples/douglas_2024/basecomputeaug.md
index 1c6114f..235e2d7 100644
--- a/examples/douglas_2024/basecomputeaug.md
+++ b/examples/douglas_2024/basecomputeaug.md
@@ -151,8 +151,9 @@ sym = 0;
 real[int] ik(sym.n), ik2(sym.n), ik3(sym.n);
 real iomega = 0.0, iomega2 = 0.0, iomega3 = 0.0;
 include "eqns.idp" // load equations
-Mat Jp(J.n, mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
-Ja = [[J,Jp],[Jp',0]]; // make dummy Jacobian
+real[int] Jp, vaug = vJp(0, XMh, tgv = -10);
+ChangeNumbering(J, vaug, Jp);
+Ja = [[J, Jp], [Jp', 0]];
 // Function to build residual operator in PETSc numbering
 func real[int] funcR(real[int]& qPETSc) {
     ChangeNumbering(J, ub[], qPETSc(0:qPETSc.n - (mpirank == 0 ? 2 : 1)), inverse = true, exchange = true);
@@ -160,9 +161,7 @@ func real[int] funcR(real[int]& qPETSc) {
     broadcast(processor(0), c);
     real[int] RPETSc, R = vR(0, XMh, tgv = TGV);
     ChangeNumbering(J, R, RPETSc);
-    ub[] .*= J.D;
-    real pavg, pavgl = int2d(Th)( y*ubp );
-    mpiAllReduce(pavgl, pavg, mpiCommWorld, mpiSUM);
+    real pavg = J(vaug, ub[]);
     if(mpirank == 0) {
         RPETSc.resize(RPETSc.n+1); // Append 0 to residual vector on proc 0
         RPETSc(RPETSc.n-1) = pavg;
@@ -175,10 +174,6 @@ func int funcJ(real[int]& qPETSc) {
     if(mpirank == 0) c = qPETSc(qPETSc.n-1);
     broadcast(processor(0), c);
     J = vJ(XMh, XMh, tgv = TGV);
-    real[int] vP, vaug = vJp(0, XMh, tgv = -10);
-    ChangeNumbering(J, vaug, vP); // FreeFEM to PETSc
-    matrix tempPms = [[vP]]; // dense array to sparse matrix
-    ChangeOperator(Jp, tempPms, parent = Ja); // send to Mat
     return 0;
 }
 // set up Mat parameters
diff --git a/examples/douglas_2024/modecomputeaug.md b/examples/douglas_2024/modecomputeaug.md
index 7caeea1..f19e25d 100644
--- a/examples/douglas_2024/modecomputeaug.md
+++ b/examples/douglas_2024/modecomputeaug.md
@@ -122,7 +122,6 @@ if (filein == "" && meshin != "") {
   ub[] = um[].re;
   filein = fileout + ".base";
 }
-Mat<complex> Jp(J.n, mpirank == 0 ? 1 : 0), ZZ(mpirank == 0 ? 1 : 0, mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
 
 sym = parsesymstr(symstr);
 complex[int] ik(sym.n), ik2(sym.n), ik3(sym.n);
@@ -144,14 +143,8 @@ complex shiftf = string2complex(targetf);
   M = vM(XMh, XMh, tgv = -20);
   complex[int] vP, vaug = vJp(0, XMh, tgv = -20);
   ChangeNumbering(J, vaug, vP); // FreeFEM to PETSc
-  matrix<complex> tempMx = [[vP]]; // dense array to sparse matrix
-  ChangeOperator(Jp, tempMx); // send to Mat
-  tempMx = [[0]];
-  ChangeOperator(ZZ, tempMx); // send to Mat
-  Mat<complex> Jatemp = [[J, Jp], [Jp', 0]]; // make dummy Jacobian
-  Mat<complex> Matemp = [[M, 0], [0, ZZ]]; // make dummy Jacobian
-  MatConvert(Jatemp, Ja);
-  MatConvert(Matemp, Ma);
+  Ja = [[J, vP], [vP', 0]];
+  Ma = [[M, 0], [0, 0]];
 }
 for (int n = 0; n < ntarget; ++n){
   if (ntarget > 1) shift = shifts + (shiftf - shifts)*real(n)/real(ntarget - 1);