Equations {D[R, y1] == 0, D[R, y2] == 0, D[R, y3] == 0}
mean that R
is a constant. Therefore, one equation needs to be solved R == R0 = const
. Consider a potential field {c1,c2,c3}=Grad[f[y1,y2,y3],{y1,y2,y3}]
, then the equation is reduced to the form
$(\vec y-\nabla f)^2+2\nabla ^2f=R0$, $\vec y=(y1,y2,y3)$
This nonlinear equation in 3D cannot be solved even numerically with the
help of NDSolve of Mathematica 11.3.
Since Mathematica 11.3 implements the finite element method for solving linear problems for elliptic equations in 3D, we use this method together with the fixed point method for solving a nonlinear problem:
R0=6; F[0][y1_, y2_, y3_] := (y1^2 + y2^2 + y3^2)/2
Do[F[i] = NDSolveValue[{(y3 -
\!\(\*SuperscriptBox[\(F[\(-1\) + i]\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3])^2 + (y2 -
\!\(\*SuperscriptBox[\(F[\(-1\) + i]\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3])^2 + (y1 -
\!\(\*SuperscriptBox[\(F[\(-1\) + i]\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3])^2 + 2 (
\!\(\*SuperscriptBox[\(f\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "2"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3] +
\!\(\*SuperscriptBox[\(f\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "2", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3] +
\!\(\*SuperscriptBox[\(f\),
TagBox[
RowBox[{"(",
RowBox[{"2", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[y1, y2, y3]) == R0,
DirichletCondition[f[y1, y2, y3] == 0, True]},
f, {y1, 0, 1}, {y2, 0, 1}, {y3, 0, 1}], {i, 1, k}];
The solution of the equation at the last step and the difference of solutions at the last two steps
{Plot3D[F[k][0.5, y2, y3], {y2, 0, 1}, {y3, 0, 1}, Mesh -> None,
ColorFunction -> Hue], Plot[F[k][.5, .5, y3] - F[k - 1][.5, .5, y3], {y3, 0, 1}]}

The distribution of the vector field {c1,c2,c3}
in space, on the plane and the line level of one component c1[y1,y2,.5]
{VectorPlot3D[
Evaluate[{D[F[k][y1, y2, y3], y1], D[F[k][y1, y2, y3], y2],
D[F[k][y1, y2, y3], y3]}], {y1, 0, 1}, {y2, 0, 1}, {y3, 0, 1},
VectorColorFunction -> Hue],
VectorPlot[
Evaluate[{D[F[k][y1, y2, y3], y1], D[F[k][y1, y2, y3], y2]} /.
y3 -> .5], {y1, 0, 1}, {y2, 0, 1}, VectorColorFunction -> Hue],
ContourPlot[
Evaluate[D[F[k][y1, y2, y3], y1] /. y3 -> .5], {y1, 0, 1}, {y2, 0,
1}, Contours -> 20, ColorFunction -> Hue]}

NDSolve
won't solve for general solutions. Add exactly enough initial conditions for your system to specify a specific solution and see what happens. $\endgroup$