# NDSolve second-order PDE with 4 variables

I'm trying to numerically solve a second-order partial differential equation in 4 dimensions. The PDE is a Fokker-Planck equation and the variables are the real and imaginary parts of the field of two coupled oscillators. I try to solve this with NDSolve using Dirichlet boundaries, but I get a strange error. I first create a 4-dimensional region in which the solution has to be found.

xmin = 0;
xmax = 13;
ymin = -10;
ymax = 10;
dbc = {DirichletCondition[P[x1, y1, x2, y2] == 0.0, True]};
Ω4d =
Parallelepiped[{xmin, ymin, xmin,
ymin}, {{xmax - xmin, 0, 0, 0}, {0, ymax - ymin, 0, 0}, {0, 0,
xmax - xmin, 0}, {0, 0, 0, ymax - ymin}}]


Then I create a very nasty PDE which I call pde $$\left(-0.02 \text{x1}^2 \text{y1}+1. \text{x1}^3+1. \text{x1} \text{y1}^2+1. \text{x1}+0.1 \text{x2}-0.02 \text{y1}^3-0.002 \text{y2}\right) P^{(0,1,0,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+\left(-1. \text{x1}^2 \text{y1}-0.02 \text{x1}^3-0.02 \text{x1} \text{y1}^2-0.002 \text{x2}-1. \text{y1}^3-1. \text{y1}-0.1 \text{y2}\right) P^{(1,0,0,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+\left(0.1 \text{x1}-0.02 \text{x2}^2 \text{y2}+1. \text{x2}^3+1. \text{x2} \text{y2}^2+1. \text{x2}+0.018 \text{y1}-0.02 \text{y2}^3-0.02 \text{y2}\right) P^{(0,0,0,1)}(\text{x1},\text{y1},\text{x2},\text{y2})+\left(0.018 \text{x1}-1. \text{x2}^2 \text{y2}-0.02 \text{x2}^3-0.02 \text{x2} \text{y2}^2-0.02 \text{x2}-0.1 \text{y1}-1. \text{y2}^3-1. \text{y2}\right) P^{(0,0,1,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+\left(-0.04 \text{x1}^2-0.04 \text{x2}^2-0.04 \text{y1}^2-0.04 \text{y2}^2-0.04\right) P(\text{x1},\text{y1},\text{x2},\text{y2})=0.02 P^{(0,0,0,2)}(\text{x1},\text{y1},\text{x2},\text{y2})+0.02 P^{(0,0,2,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+0.002 P^{(0,2,0,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+0.002 P^{(2,0,0,0)}(\text{x1},\text{y1},\text{x2},\text{y2})$$ I try to solve this with

NDSolve[{pde, dbc}, P, {x1, y1, x2, y2} ∈ Ω4d]


But I get the following error message

NDSolve::fem3: -- Message text not found -- (Parallelepiped[{0,-10,0,-10},{{13,0,0,0},{0,20,0,0},{0,0,13,0},{0,0,0,20}}])


If I solve with

NDSolve[{pde, dbc}, P, {x1, xmin, xmax}, {y1, ymin, ymax}, {x2, xmin,xmax}, {y2, ymin, ymax}]


I get

NDSolve::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.


Does anyone know what I'm doing wrong here?