Return to Revisions

Here's my FDM-based solution for your problem:

{{xl, xr}, {yl, yr}} = {{0, 33}, {0, 15}};
sf = 2;
nx = sf xr; ny = sf yr;
dx = (xr - xl)/nx; dy = (yr - yl)/ny;
xmidl = 15; xmidr = 18; ymid = 9;
h1 = 14; h2 = 2;
formula = Select[
   Flatten@Table[
      h[x - dx, y] + h[x + dx, y] + h[x, y - dy] + h[x, y + dy] - 4 h[x, y] == 0, {x, xl,
        xr, dx}, {y, yl, yr, dy}][[2 ;; -2, 2 ;; -2]], 
   FreeQ[#, h[x_, y_] /; xmidl < x < xmidr && y > ymid] &];
oneSideD1[most__, "left"] := oneSideD1[most, -1]
oneSideD1[most__, "right"] := oneSideD1[most, 1]
oneSideD1[h_, x_, step_, direction : 1 | -1] := 
 direction ((3 h@x)/(2 step) - (2 h[x - direction step])/step + h[x - 2 direction step]/
    step)
bcxl = Table[oneSideD1[h[#, y] &, xl, dx, "left"] == 0, {y, yl, yr, dy}][[2 ;; -2]];
bcxr = Table[oneSideD1[h[#, y] &, xr, dx, "right"] == 0, {y, yl, yr, dy}][[2 ;; -2]];
bcyl = Table[oneSideD1[h[x, #] &, yl, dy, "left"] == 0, {x, xl, xr, dx}];
bcyr@1 = Table[h[x, yr] == h1, {x, xl, xmidl, dx}];
bcyr@5 = Table[h[x, yr] == h2, {x, xmidr, xr, dx}];
bcyr@3 = Table[
    oneSideD1[h[x, #] &, ymid, dy, "right"] == 0, {x, xmidl, xmidr, dx}][[2 ;; -2]];
bcyr@2 = Table[
    oneSideD1[h[#, y] &, xmidl, dx, "right"] == 0, {y, ymid, yr, dy}][[2 ;; -2]];
bcyr@4 = Table[
    oneSideD1[h[#, y] &, xmidr, dx, "left"] == 0, {y, ymid, yr, dy}][[2 ;; -2]];
set = Flatten@{formula, bcxl, bcxr, bcyl, bcyr /@ Range@5};
var = Union@Cases[set, h[a_, b_], ∞];
{b, mat} = CoefficientArrays[set, var];
sol = LinearSolve[mat, -N@b];
coord = List @@@ var;
ListPointPlot3D[Flatten /@ ({coord, sol}\[Transpose]), PlotRange -> All]

Remark

1. I've used one-sided difference formula i.e. oneSideD1 to discretize the Neumann boundary condtions. For your simple Neumann b.c., it's not a bad idea to handle them with reflection of course, but do notice reflection is hard to extend to more general cases.

2. sf should be an Integer.

3. You can also use

    SetAttributes[h, NHoldAll]    
    sol2 = Solve[N@set, var]; // AbsoluteTiming
   

to solve the equation set, but this approach is slower. (The speed difference isn't obvious in this simple case though.)