I have developed the code below to solve two PDEs; first mu[x,y] is solved for, then the results of mu are used to solve for phi[x,y]. The code works and converges on a solution as is, however, I would like to decrease the size of a, b, and d even further. To accurately represent the physical process I am trying to simulate, a, b, and d would need to be ~100-1000x smaller. If I make them smaller, I don’t believe the solution has actually converged because the values for phi along the right boundary change significantly with a change in mesh size (i.e. if I make them smaller and the code below produces a value of phi=-0.764 at the midpoint between y2 and y3 along the right boundary, a change in size1 to 10^-17 and size2 to 10^-15, changes that value of phi to -0.763, and a change in size2 to 10^-16 changes that value again to -0.860), but I cannot make the mesh size any smaller without Mathematica crashing.

Are there any better ways to create the mesh that would be less computationally taxing and allow it to be more refined in the regions of interest? Or are there any ways to make the code in general less computationally expensive so that I can further refine the mesh?

`ClearAll["Global`*"] Needs["NDSolve`FEM`"] (* 1) Define Constants*) e = 1.60217662*10^-19; F = 96485; kb = 1.381*10^-23; sigi = 18; sigini = 0; sigeni = 2*10^6; T = 1000; n = -0.02; c = 1; pH2 = 0.2; pH2O = 1 - pH2; pO2 = 1.52*^-19; l = 10*10^-6; a = 100*10^-7; b = 50*10^-7; d = 300*10^-7; y1 = 0.01; y2 = 0.5*y1; y3 = y2 + a; y4 = y3 + d; y5 = y4 + b; mu1 = 0; mu2 = -5.98392*^-19; phi1 = 0; (* 2) Create mesh*) m = 0.1*l; size1 = 10^-16; size2 = 10^-15; size3 = 10^-7; mrf = With[{rmf = RegionMember[ Region@RegionUnion[Disk[{l, y2}, m], Disk[{l, y3}, m], Disk[{l, y4}, m], Disk[{l, y5}, m]]]}, Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices]; Which[rmf[{x, y}], area > size1, (0 <= x <= l && y2 - l <= y <= y2 + l), area > size2, (0 <= x <= l && y3 - l <= y <= y3 + l), area > size2, (0 <= x <= l && y4 - l <= y <= y4 + l), area > size2, (0 <= x <= l && y5 - l <= y <= y5 + l), area > size2, True, area > size3]]]]; mesh = DiscretizeRegion[Rectangle[{0, 0}, {l, y1}], MeshRefinementFunction -> mrf]; (* 3) Solve for mu*) bcmu = {DirichletCondition[mu[x, y] == mu1, (x == 0 && 0 < y < y1)], DirichletCondition[ mu[x, y] == mu2, (x == l && y2 <= y <= y3) || (x == l && y4 <= y <= y5)]}; solmu = NDSolve[{Laplacian[mu[x, y], {x, y}] == 0 + NeumannValue[0, y == 0 || y == y1 || (x == l && 0 <= y < y2) || (x == l && y3 < y < y4) || (x == l && y5 < y < y1)], bcmu}, mu, {x, y} \[Element] mesh, WorkingPrecision -> 50]; (* 4) Solve for electronic conductivity everywhere*) pO2data = Exp[(mu[x, y] /. solmu)/kb/T]; sige0 = 2.77*10^-7; sigedata = Piecewise[{{sige0*pO2data^(-1/4), 0 <= x <= l - m}, {sige0*pO2data^(-1/4), (l - m < x <= l && 0 <= y < y2)}, {(sigeni - sige0*(pO2data /. x -> l - m)^(-1/4))/m*(x - (l - m)) + sige0*(pO2data /. x -> l - m)^(-1/4), (l - m < x <= l && y2 <= y <= y3)}, {sige0*pO2data^(-1/4), (l - m < x <= l && y3 < y < y4)}, {(sigeni - sige0*(pO2data /. x -> l - m)^(-1/4))/m*(x - (l - m)) + sige0*(pO2data /. x -> l - m)^(-1/4), (l - m < x <= l && y4 <= y <= y5)}, {sige0*pO2data^(-1/4), (l - m < x <= l && y5 < y <= y1)}}]; (* 5) Solve for phi*) Irxn = -(2*F)*(c*pO2^n ); A = (Irxn - sigi/(4*e)*(D[mu[x, y] /. solmu, x] /. x -> l))/(-sigi); B = sigi/(4*e)*(D[mu[x, y] /. solmu, x] /. x -> l)/(sigi + sigedata /. x -> l - m); bcphi = DirichletCondition[phi[x, y] == phi1, (x == 0 && 0 < y < y1)]; solphi = NDSolve[{Laplacian[phi[x, y], {x, y}] == 0 + NeumannValue[0, y == 0 || y == y1 || (x == l && 0 <= y < y2) || (x == l && y3 < y < y4) || (x == l && y5 < y < y1)] + NeumannValue[-A[[1]], (x == l && y2 <= y <= y3)] + NeumannValue[-B[[1]], (x == l && y4 <= y <= y5)], bcphi}, phi, {x, y} \[Element] mesh, WorkingPrecision -> 50]; (* 6) Print values to check for convergence*) P[x_, y_] := phi[x, y] /. solphi; P[l, (y3 - y2)/2 + y2] P[l, (y5 - y4)/2 + y4] `