In the case of two metal objects, we can set the potential of each object as $V_1, V_2$. Then the code for a numerical solution in 2D is Needs["NDSolve`FEM`"]; (*Define Boundaries*)air = Rectangle[{-5, -5}, {5, 5}]; object1 = Disk[]; object2 = Rectangle[{2, 0}, {2.5, 2}]; reg12 = RegionUnion[object1, object2]; reg = RegionDifference[air, reg12]; mesh = ToElementMesh[reg, MeshRefinementFunction -> Function[{vertices, area}, area > 0.001 (0.1 + 10 Norm[Mean[vertices]])]] mesh["Wireframe"] eq = Laplacian[u[x, y], {x, y}]; V1 = 1; V2 = -2; bc = {DirichletCondition[u[x, y] == V1, x^2 + y^2 == 1], DirichletCondition[ u[x, y] == V2, (x == 2 || x == 2.5 && 0 <= y <= 2) || (y == 0 || y == 2 && 2 <= x <= 2.5)]}; U = NDSolveValue[{eq == 0, bc}, u, {x, y} ∈ mesh]; ef = -Grad[U[x, y], {x, y}]; Visualisation of solution {DensityPlot[U[x, y], {x, y} ∈ reg, ColorFunction -> "Rainbow", PlotLegends -> Automatic, FrameLabel -> Automatic, PlotPoints -> 50, PlotRange -> {{-4, 4}, {-4, 4}}], StreamDensityPlot[Evaluate[ef], {x, y} ∈ reg, ColorFunction -> "Rainbow", PlotLegends -> Automatic, FrameLabel -> {x, y}, StreamStyle -> LightGray, VectorPoints -> Fine, PlotRange -> {{-1, 3}, {-1, 3}}]} [![Figure 1][1]][1] Update 1. Next code is devoted to solve electrostatic problem for combination of dielectric and conducting objects (glass cylinder and metal strip). For dielectric we put electric charge $q_1$, and for metal we put potential $V_2$. Code: Needs["NDSolve`FEM`"]; par = {eps1 -> 3.5, eps2 -> 1.0}; air = Rectangle[{-5, -5}, {5, 5}]; object1 = Disk[]; q1 = 1; vol1 = NIntegrate[1, {x, y} ∈ object1]; rho1 = q1/vol1; object2 = Rectangle[{2, 0}, {2.5, 2}]; rho[x_, y_] := rho1 Boole[{x, y} ∈ object1]; eps[x_, y_] := eps2 + (eps1 - eps2) Boole[{x, y} ∈ object1]; reg = RegionDifference[air, object2]; mesh = ToElementMesh[reg, MeshRefinementFunction -> Function[{vertices, area}, area > 0.001 (0.1 + 10 Norm[Mean[vertices]])]] mesh["Wireframe"] V2 = -2; eq = Inactive[Div][ eps[x, y] Inactive[Grad][u[x, y], {x, y}], {x, y}] == -2 Pi rho[x, y]; bc = DirichletCondition[u[x, y] == V2, {x, y} ∈ object2]; U = NDSolveValue[{eq /. par, bc}, u, {x, y} ∈ mesh]; ef = -Grad[U[x, y], {x, y}]; Visualisation {DensityPlot[U[x, y], {x, y} ∈ mesh, ColorFunction -> "Rainbow", PlotLegends -> Automatic, FrameLabel -> Automatic, PlotPoints -> 50, PlotRange -> {{-4, 4}, {-4, 4}}], StreamDensityPlot[Evaluate[ef], {x, y} ∈ reg, ColorFunction -> "Rainbow", PlotLegends -> Automatic, FrameLabel -> {x, y}, StreamStyle -> LightGray, VectorPoints -> Fine, PlotRange -> {{-1, 3}, {-1, 3}}], StreamDensityPlot[Evaluate[ef], {x, y} ∈ reg, ColorFunction -> Hue, FrameLabel -> {x, y}, StreamStyle -> Blue, PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}}, PlotLegends -> Automatic]} [![Figure 2][2]][2] [1]: https://i.sstatic.net/uY6Ey.png [2]: https://i.sstatic.net/Aavj7.png