# Solving an 'odd' differential equation with NDSolve

I need to solve a differential equation of the type

$\qquad \partial_{x_1}y(x_1,x_2)= y(x_1,x_1)\,y(x_1,x_2)$

with initial condition

$\qquad y(0,x_2)=x_2$.

Now if I try to solve this with NDSolve I get an error that tells me "the arguments should be ordered consistently" (NDSolve::conarg), here is my code:

s =
NDSolve[
{D[y[x1, x2], x1] == y[x1, x1] y[x1, x2], y[0, x2] == x2},
y, {x1, 0, 1}, {x2, 0, 1}]


On my daylong search I haven't found any helpful solution for this seemingly easy problem.

• Since it's not mentioned in the question, I'll point out that the cause of the error is that the arguments of y[x1, x1] and y[x1, x2] are not the same. Commented Sep 10, 2018 at 16:35
• The y[x1,x1] is the point of the problem, otherwise it would be trivial. Sorry, I should have clarified that. I want it to work this way. Commented Sep 10, 2018 at 17:16
• Just wanted to make sure the meaning of the error was clear. I wonder, what makes this problem seem easy? Is it a standard type? Are there standard methods for numerically integrating it? Commented Sep 10, 2018 at 22:09
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Commented Sep 12, 2018 at 3:55

This PDE can be solved symbolically. To begin, replace y[x1, x1] by f[x1], which DSolve can handle without difficulty.

DSolve[{D[y[x1, x2], x1] == f[x1] y[x1, x2], y[0, x2] == x2},
y[x1, x2], {x1, x2}] [[1, 1]]
(* y[x1, x2] -> E^(Integrate[-f[K[1]], {K[1], 1, 0}]
- Integrate[-f[K[1]], {K[1], 1, x1}])*x2


From this solution, it is clear that y[x1, x2]/x2 is independent of x2. Consequently, y[x1, x1]/x1 == y[x1, x2]/x2. Substituting this identity into the PDE then yields,

s = DSolve[{D[y[x1, x2], x1] == (x1/x2) y[x1, x2]^2, y[0, x2] == x2},
y[x1, x2], {x1, x2}] [[1, 1]]
(* y[x1, x2] -> -((2 x2)/(-2 + x1^2)) *)

Plot3D[s//Last, {x1, 0, 1}, {x2, 0, 1}, AxesLabel -> {x1, x2, y},
ImageSize -> Large, LabelStyle -> {Bold, Black, Medium}]


Not surprisingly, the plot closely resembles the numerical result of Michael E2. Note that this method works with this PDE for any initial condition.

This solution can, if desired, be validated by back-substitution:

{D[y[x1, x2], x1] == y[x1, x1] y[x1, x2], y[0, x2] == x2} /.
y -> Function[{x1, x2}, -((2 x2)/(-2 + x1^2))]
(* {True, True} *)

• Writing down f[x1], you implicitly claim it a known function, but it is not so. Therefore, your solution is built on the sand. Am I not right? Commented Sep 12, 2018 at 4:57
• @user64494 No. I do not claim that f[x1] is known, only that y[x1, x1] is a function only of x1. Commented Sep 12, 2018 at 5:02
• OK, Function[{x1, x2}, -((2 x2)/(-2 + x1^2))] is a solution. I still have doubts concerning your DSolving. Commented Sep 12, 2018 at 5:20

One could try a iterative approach approximating y[x1, x1] with the previous, if the iterations converge to solution:

foo::div = "Error increased. Divergent?";
foo::slwcon = "Slow convergence. Error  after  steps.";
tol = 10^-6;
NestWhile[
Function[{s1, s2, err1, err2, iter},
With[{s =
NDSolveValue[{D[y[x1, x2], x1] == s1[x1, x1] y[x1, x2],
y[0, x2] == x2}, y, {x1, 0, 1}, {x2, 0, 1}]},
{s, s1,
Flatten[s1["ValuesOnGrid"], 1] - s @@@ Flatten[s1["Grid"], 1] //
Abs // Max,
err1,
iter + 1}
]] @@ # &,
{NDSolveValue[{D[y[x1, x2], x1] == 1*y[x1, x2], y[0, x2] == x2},
y, {x1, 0, 1}, {x2, 0, 1}],
1 &,
10., 100., 0},
Function[{s1, s2, err1, err2, iter},
If[err1 > err2,
Message[foo::div]; False,
If[err1 > tol,
If[iter == 10, Message[foo::slwcon, err1, iter]];
True,
False]
]
] @@ # &,
1, 100
]
s = {y -> First[%]};


Plot3D[y[x1, x2] /. s, {x1, 0, 1}, {x2, 0, 1}, AxesLabel -> Automatic]


The error is so-so:

Plot3D[D[y[x1, x2], x1] - y[x1, x1] y[x1, x2] /. s // Evaluate,
{x1, 0, 1}, {x2, 0, 1}, AxesLabel -> Automatic]


• Thanks for your effort! Commented Sep 11, 2018 at 14:38