I have a solution for the Laplace equation with simple Dirichlet boundary conditions. The top side of the square is kept constant at 100 while other sides are fixed at 0. The code is here:
ClearAll["Global`*"] Needs["NDSolve`FEM`"] square = Rectangle[{0, 0}, {10, 10}]; mesh = ToElementMesh[DiscretizeRegion[square], MaxCellMeasure -> 0.01]; bc = { DirichletCondition[u[x, y] == 0, y == 0 || x == 10 || x == 0], DirichletCondition[u[x, y] == 100, y == 10] }; sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0 , bc}, u, {x, y} \[Element] mesh]; DensityPlot[sol[x, y], {x, y} \[Element] mesh, Mesh -> None, ColorFunction -> "Rainbow", PlotRange -> All, PlotLegends -> Automatic] Here is the plotted solution:
Now what I would like to do is to set Neumann boudary conditions for the left and right walls. I want to fix the gradient at 0 for these two walls.
I approach it by getting rid of || x == 10 || x == 0 in the old boundary condition setup and introducing:
Derivative[1,0][u][0,y] == 0, Derivative[1,0][u][10,y] == 0 This doesn't work for me. Could you please help me find the mistake?


|| x == 10 || x == 0and it should be OK. $\endgroup$