# Problems using NDSolveValue to solve system of PDEs over 3D region

I am trying to solve a system of partial differential equations over a cube region. The text is as follows:

<< NDSolveFEM
\[Beta] =
Parallelepiped[{-1, -1, -1}, {{0, 0, 2}, {2, 0, 0}, {0, 2, 0}}]
mesh = ToBoundaryMesh[\[Beta]]
mesh2 := ToElementMesh[mesh, MaxCellMeasure -> 0.01]
J = 10;
w = 1;
Y = 10;
{Asol, Bsol, psol} =
NDSolveValue[{D[A[x, y, z], z] - D[p[x, y, z], x] == 0,
D[p[x, y, z], y] - D[B[x, y, z], z] == 0,
D[B[x, y, z], x] - D[A[x, y, z], y] == 1,
DirichletCondition[p[x, y, z] == 0, z == -1],
DirichletCondition[p[x, y, z] == 0, z == 1],
DirichletCondition[B[x, y, z] == 0, y == -1],
DirichletCondition[A[x, y, z] == 0, x == -1],
DirichletCondition[A[x, y, z] == 0, x == 1],
DirichletCondition[B[x, y, z] == 0, y == 1]}, {A, B,
p}, {x, y, z} \[Element] mesh2]


Everything seems to work fine, however, when I evaluate the value of D[B[x, y, z], x] - D[A[x, y, z], y] at any point, I am returned with 0, which does not match the equation that I entered, that estipulates that D[B[x, y, z], x] - D[A[x, y, z], y] == 1. So essentialy, mathematica is returning me a solution to my system of PDEs that doesn't actually solve the system. How can I fix this ?

• The solution of the system of equations diverges. The task is incorrectly set. – Alex Trounev Nov 13 '18 at 22:38