# Using NDSolve on piecewise-defined system of ODEs

I'd like to use NDSolve to solve a system of n coupled 1st-order linear ODEs expressed in matrix-vector form, but which are definied piecewise in n-space (in particular, to represent different linear regimes of a piecewise-linear state-space model). I can successfully use NDSolve in a case such as

sol = X[t] /. NDSolve[{
X'[t] == A . X[t],
X[0] == X0}, X, {t, 0, 10}];


where A is some matrix is X0 is the initial condition. However, I'd like different behavior at different areas in space per a piecewise function, namely

sol1 = X[t] /. NDSolve[{
X'[t] == xdot[X[t]],
X[0] == X0}, X, {t, 0, 10}];


where, for example

xdot[X_] := Piecewise[{{A2 .
X, ((X[[1]] - X[[3]] > 0) && (X[[3]] >=
d/2)) || ((X[[1]] - X[[3]] < 0) && (X[[3]] <= -d/2))}}, A1 . X]


but this fails because the piecewise function must use part specifications but they fail when applied to the undefined 'X[t]' used in NDSolve.

Looking for ways to remedy this or suggestions for other ways to go about implementing the model. Thanks!

• Have you tried Indexed[X, 3]? Mar 24, 2023 at 23:54
• @MichaelE2 ah, that worked, thanks! Got a solution - didn't behave quite as expected, though I'd probably open a separate issue/question for that. Mar 25, 2023 at 5:42
• What is A and X0? Please post the complete code. Mar 25, 2023 at 10:42

Its always a good idea to introduce all variables in a transparent, adjustable way. In order to define functions with conditions never use Piecewise. Make general and conditional definitions.

Consider the oscillator equation with a breaking coil

Force[x_] := -x /; Abs[x] < 1
Force[x_] := 0 /; Abs[x] > 1

Plot[Force[x], {x, -4, 4}]

h[t_] = g[t] /.
NDSolve[Evaluate[{g''[t] == Force[g[t]], g[0] == 0,
g'[0] == 1.002}], g[t], {t, -3, 3}]

Plot[h[t], {t, -4, 4}]


In your case try (I didn't)

xdot[A_, B_, X_, d_] := B . X
xdot[A_, B_, X_, d_] :=
A . X /; ((X[[1]] - X[[3]] > 0) && (X[[3]] >=
d/2)) || ((X[[1]] - X[[3]] < 0) && (X[[3]] <= -d/2))

A1 = RandomReal[{-1, 1}, {2, 6}];

A2 = RandomReal[{-1, 1}, {2, 6}];

X = RandomReal[{-12, 12}, {6}];
Plot[Total[xdot[A1, A2, X, d]], {d, -120, 120}]