**it** can read

not much space left. Take an unit square instead for the independent variables.

One gets explanations for the messages in the help ref/message/NDSolve/deqn, so a derivative on the boundary is not specified by

D[z[x,0],y]==0

because that gives True (z[x,0] does not depend on y), but using

(D[z[x,y],y]/.y->0)==0

if done so one has

In[108]:= NDSolveValue[{D[z[x, y], {x, 4}] +

2 D[z[x, y], {x, 2}, {y, 2}] + D[z[x, y], {y, 4}] == q/d,

z[x, 0] == 0,

z[0, y] == 0, (D[z[x, y], y] /. y -> 0) ==

0, (D[z[x, y], x] /. x -> 0) == 0 }, z, {x, 0, 1}, {y, 0, 1}]

During evaluation of In[108]:= NDSolveValue::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable. >>

now one reads ref/message/NDSolve/ivone and considers the type of the equation, especially it states

The partial differential equation solver used by NDSolve is designed to solve initial value problems.

purely elliptic problems are not in the scope of NDSolve[], if I understand it right.