# Density plot of asymptotic behaviour of a linear system

Good morning,

I have an ODE linear system $\dot{\vec f} = A \vec f$ depending parametrically to two real parameters, s and p. I'd like to get a density map of the solutions of $f_A/(f_A+f_B)$ at $t=20$, in the space of parameters ${s, p} \in [0; 1]$.

I post my try:

rA = 2.1
rB = 1.0
fA := ParametricNDSolveValue[{A'[ t] == (rA - 1 + p) A[ t] + s B[t], B'[t] == (1 - p) A[t] + (rB - s) B[t], A[0] == 2, B[0] == 2}, A, {t, 0, 20}, {{s, p} \[Element] Reals}]
fB := ParametricNDSolveValue[{A'[ t] == (rA - 1 + p) A[ t] + s B[t], B'[t] == (1 - p) A[t] + (rB - s) B[t], A[0] == 2, B[0] == 2}, B, {t, 0, 20}, {{s, p} \[Element] Reals}]
DensityPlot[fA[s][p]/(fA[s][p] + fB[s][p]), {s, 0, 1}, {p, 0, 1}]


The plot I obtain can't be what I need, as I never ask Mathematica to evaluate the functions at $t=20$. I'm really new to Mathematica and I've the impression to be missing something really basic about its semantics. I hope to learn from your help, thank you in advance!

-
I know that nowhere in my code there is $t=20$, but when I tried few things to do that I got errors like "ParametricNDSolve::dsvar: 20 cannot be used as a variable" –  Orso Jul 8 '14 at 9:06
I cut and pasted your code and the plot works? Try Clear[t] first, it looks like you may have executed t=20 before. –  Ymareth Jul 8 '14 at 9:16
Hi @Ymareth and thank you for your help. Indeed, I already do Clear[t] before the lines I posted but it doesn't change anything in my case... –  Orso Jul 8 '14 at 9:21
You're right, I pasted it into a new notebook and I obtain a plot. The problem is that nobody told the machine to plot at $t=20$, so I guess the plot is not at all what I need... –  Orso Jul 8 '14 at 9:38

I may have misunderstood the aim. I post this as motivating:

I have re-written the system of equations (no guarantee error free):

mat = {{1.1 + p, s}, {1 - p, 1 - s}}; pf =
ParametricNDSolveValue[{{a'[t], b'[t]} ==
mat.{a[t], b[t]}, {a[0], b[0]} == {2, 2}}, {a[t], b[t]}, {t, 0,
20}, {s, p}];
foi[par1_, par2_, u_] := (#1/(#2 + #1) & @@ pf[par1, par2]) /. t -> u


Looking at the solutions $a(t),b(t)$:

Manipulate[
Column[{ParametricPlot[pf[v1, v2] /. t -> w, {w, 0, 20},
PlotRange -> {{0, 10^12}, {0, 10^12}}, Frame -> True,
FrameLabel -> {"a", "b"}],
LogPlot[Evaluate[pf[v1, v2] /. t -> w], {w, 0, 20},
PlotLegends -> {"a(t)", "b(t)"}
, PlotRange -> {0, 10^12}, Frame -> True,
FrameLabel -> {"t", "a|b"}]}], {{v1, 0, "s"}, 0,
1}, {{v2, 0, "p"}, 0, 1}]


As expected with these linear system with constant positive coefficients->exponentials

The desired density plot for function (foi evaluated at t=20) above:

DensityPlot[foi[x, y, 20], {x, 0, 1}, {y, 0, 1},
ColorFunction -> "SunsetColors", PlotLegends -> Automatic,
MeshFunctions -> (#3 &), Mesh -> {Range[0, 1, 0.1]},
FrameLabel -> {"s", "p"}]


or in 3D:

Plot3D[foi[x, y, 20], {x, 0, 1}, {y, 0, 1},
ColorFunction -> (ColorData["SunsetColors"][#3] &),
MeshFunctions -> (#3 &), Mesh -> {Range[0, 1, 0.1]},
AxesLabel -> {"s", "p",
"\!$$\*FractionBox[\(a \((20)$$\), $$a \((20)$$ + b \
$$(20)$$\)]\)"}]


Apologies if I have misinterpreted your aim.

-
Exactly, @ubpqdn, that's what I wanted! Thank you so much, I've learnt a lot from your solution! –  Orso Jul 8 '14 at 13:48