# Square root problem

Consider the following code

NDSolveValue[{
x'[t] == x[t]^a,
x[0] == 0
}, x, {t, 0, 10}]


which solves the equation $$x'(t)=x(t)^a$$, for some exponent $$a$$, with $$x(0)=0$$. If I want to solve two identical equations, I can simply add a new initial condition

NDSolveValue[{
x'[t] == x[t]^a,
x[0] == {0, 0}
}, x, {t, 0, 10}]


I want to solve this system when $$a=1/2$$ (square root), however this doesn't work for more than one equation and I get, among others, the error

"NDSolveValue::nbnum1: The function value {0,0}==0 is not True or False when the arguments are {0.,{0.,0.},{0.,0.}}."


Why is that? If I pick any $$a\geq 1$$ it works as intended. Why is it suddenly interpreting the right hand side as a scalar for $$a<1$$? In the one equation case, square root works. Any ideas?

Edit and motivation: I want to mimic the following code

NDSolveValue[{x1'[t] == Sqrt[x1[t]], x2'[t] == Sqrt[x2[t]],
x1[0] == x2[0] == 0}, {x1, x2}, {t, 0, 10}]


and the reason I'm using this is because I would like to generate a big list of equations. For $$a\geq 1$$, at least, the main difference between this approach and mine is that here I get a list of two interpolating functions with scalar outputs, and in my case I get one interpolating function with 2 output dimensions.

In other words, what is the best way to solve $$x_1'(t)=\sqrt{x_1(t)}, x_2'(t)=\sqrt{x_2(t)}, \cdots,x_n'(t)=\sqrt{x_n(t)}\\$$ for an arbitrary $$n$$? The main reason I'm interested in this approach is because it's what I use to solve a system like this $$\begin{pmatrix} x_1'(t)\\ x_2'(t)\\ \vdots\\ x_n'(t) \end{pmatrix}=A f\begin{pmatrix} x_1(t)\\ x_2(t)\\ \vdots\\ x_n(t) \end{pmatrix}$$ where $$A$$ is some $$n\times n$$ matrix and $$f:\mathbb{R}^n\to\mathbb{R}^n$$. In this specific case, $$A$$ is the identity matrix and $$f$$ maps the square root.

• "If I want to solve two identical equations, I can simply add a new initial condition" - I don't understand that at all, as I don't understand what you mean by x[0] = {0, 0} in your initial conditions. I suspect that NDSolveValue is having the same problem interpreting what you mean. Also, exactly what code generates the error you saw? Commented Apr 19, 2021 at 14:48
• @MarcoB Here's my motivation: Imagine I want to solve a system of $n$ ODEs that are dependent on each other (say we're measuring the signal of $n$ cells interacting in a tissue, for example). I noticed that a 'compact' way of doing that in Mathematica, when the system of ODEs is really big, is to define an adjacency matrix and treating all functions as a list. For example, see what I did in this question: mathematica.stackexchange.com/questions/224470/… Also, see the edit section, to see what I mean Commented Apr 19, 2021 at 14:55
• It sounds like what you want to do is to solve the ODE for many different initial conditions "all at once"? If so, you can do something like: a = 1/2; NDSolveValue[{x'[t] == x[t]^a, x[0] == #}, x, {t, 0, 10}] & /@ Range[0, 1, 0.1] where the Range is just used here to specify all the different initial conditions. Commented Apr 19, 2021 at 14:55
• I don't understand. Even your very first system: NDSolveValue[{x'[t] == x[t]^a, x[0] == 0}, x, {t, 0, 10}] does not work without assigning a value to $a$. Commented Apr 19, 2021 at 14:58
• I think the problems is that Sqrt is a multivalued function. If you use Surd (this gives a real value) instead of Sqrt, then the following works: NDSolveValue[{x'[t] == Surd[x[t], 2], x[0] == {0, 1}}, x, {t, 0, 10}] Commented Apr 19, 2021 at 15:22

For a = 1/2, the problem is that the discontinuity processing automatically handles the singularity in the x[t]^a power (Sqrt[x[t]]) by setting up an event when x[t] == 0. Unfortunately, it doesn't notice that x[t] is a vector quantity, so that the event testing causes a slew of errors. The most telling one is NDSolveValue::ecboo, which mentions an "event" and clearly points to some trouble with the discontinuity processing.
A fix is to turn "DiscontinuityProcessing" off:
ClearAll[x, t, a];