# NDSolve::femcmsd: The spatial derivative order of the PDE may not exceed two

I was in the impression that V12 can now handle higher order PDES but for this particular case it fails to produce any solution.

pde1 = D[F[x, y], {x, 3}] + F[x, y]*D[F[x, y], {x, 2}] ==
y*(D[F[x, y], {x, 1}]*D[D[F[x, y], {x, 1}], y] -
D[F[x, y], {x, 2}]*D[F[x, y], {y, 1}]);

pde2[Pr_, L1_] =
1/Pr*D[T[x, y], {x, 2}] + F[x, y]*D[T[x, y], {x, 1}] -
2*L1*D[F[x, y], {x, 1}]*T[x, y] ==
y*(D[F[x, y], {x, 1}]*D[T[x, y], {y, 1}] -
D[T[x, y], {x, 1}]*D[F[x, y], {y, 1}]);

With[{lb = 5}, bcs = {{F[0, y] + y == -y*Derivative[0, 1][F][0, y],
Derivative[1, 0][F][0, y] == 0, T[0, y] == 1},
{Derivative[1, 0][F][lb, y] == 1, T[lb, y] == 0}}];

Clear@solfunc
With[{lb = 5},
solfunc[Pr_, L1_: 0.5] :=NDSolve[{pde1, pde2[Pr, L1], bcs}, {F, T}, {x, 0, lb}, {y, 0, 1}]]

(sollst[#] = solfunc[#]) & /@ {0.7, 3}


The above system has been solved using truncation method getting a series solution in this paper. The approach suggested by the authors is long and a tedious one. I was wondering whether it is possible to solve this system directly or there is some other efficient way to get the solution.

Any suggestion on how to solve such system?

• How did you get this system? Jan 4, 2020 at 21:57
• The system has five boundary conditions in x, as expected. But, what are the two boundary conditions in y? Jan 5, 2020 at 3:05
• With those initial conditions, it becomes an initial-value problem. Try specifying that NDSolve use the numerical method of lines. If that does not work either, descritize the PDEs in x by hand to create a large system of ODEs in y which almost certainly can be solved numerically. Jan 5, 2020 at 3:15
• Too bad. It appears you need to discretize in x by hand and then solve the resulting system of ODEs. Although tedious, it should work, I think. (The cross-derivatives in the first PDE could be an issue.) Jan 5, 2020 at 3:22
• Have a look at the message ref page ref/message/NDSolve/femcmsd. Jan 6, 2020 at 7:45

Here is a partial solution. Begin by noting that

• pde1 does not involve T and so can be solved for F and the solution then inserted into pde2 to solve for T.
• The boundary condition F[0, y] + y == -y*Derivative[0, 1][F][0, y] should be integrated to yield -y/2.
• The boundary condition Derivative[1, 0][F][lb, y] == 1 is inconsistent with the initial condition F[x,0] == 0 at x = lb, y = 0 but can be made consistent with the standard trick, Derivative[1, 0][F][lb, y] == 1 - Exp[-10 y],
• The error message cited in the question is associated with FEM and can be avoided by using Method -> "MethodOfLines".

With these changes, pde1 can be solved as follows:

lb = 5;
pde1 = D[F[x, y], {x, 3}] + F[x, y]*D[F[x, y], {x, 2}] ==
y*(D[F[x, y], {x, 1}]*D[D[F[x, y], {x, 1}], y] - D[F[x, y], {x, 2}]*D[F[x, y], {y, 1}]);
bcs = {F[0, y] == -y/2, Derivative[1, 0][F][0, y] == 0,
Derivative[1, 0][F][lb, y] == 1 - Exp[-10 y], F[x, 0] == 0};
sol = NDSolveValue[{pde1, bcs}, F[x, y], {x, 0, lb}, {y, 0, 1}, Method -> "MethodOfLines"];
Plot3D[sol, {x, 0, lb}, {y, 0, 1}, ImageSize -> Large,
AxesLabel -> {x, y, F}, LabelStyle -> {15, Black, Bold}]


• Now it is working like a charm. Thx as always for your valuable answers.
– zhk
Jan 8, 2020 at 1:08