{{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]

{{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]/
    (2 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]

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
   

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. (If you want to learn more about one-sided formula, start from page 6 of this book. )

2. sf should be an Integer.

3. You can also use

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