Define the geometry of the objects using rectangles and circles.
Combine the geometries by summing them to create a single geometry that includes all the objects.
Create a PDE model, specify coefficients, apply boundary conditions, generate a mesh, solve the PDE, and visualize the results.
I have modified the code a bit
rect1 = [3; 4; -1; 1; 1; -1; 0; 0; -0.5; -0.5];
C1 = [1; 0.7; -0.25; 0.15];
C2 = [1; -0.7; -0.25; 0.15];
C1 = [C1; zeros(length(rect1) - length(C1), 1)];
C2 = [C2; zeros(length(rect1) - length(C2), 1)];
gd = [rect1, C1, C2];
ns = char('rect1', 'C1', 'C2');
ns = ns';
sf = 'rect1 + C1 + C2';
[dl, bt] = decsg(gd, sf, ns);
model = createpde();
% Step 4: Remove overlapping regions and create a valid geometry
[dl2, bt2] = csgdel(dl, bt);
% Step 5: Generate a mesh
geometryFromEdges(model, dl2);
mesh = generateMesh(model, 'Hmax', 0.25);
% Step 6: Specify coefficients
specifyCoefficients(model, 'm', 0, 'd', 0, 'c', 1, 'a', 0, 'f', 1);
% Step 7: Apply boundary conditions
applyBoundaryCondition(model, 'dirichlet', 'Edge', 1:model.Geometry.NumEdges, 'u', 10);
% Step 8: Solve the PDE
results = solvepde(model);
% Step 9: Visualize the results
figure;
pdeplot(model, 'XYData', results.NodalSolution);
axis equal;
Do let me know if it is working for you?