Hi Otto and Pat, and thanks to both of you.

Otto; Sorry I didn't comment on (t-tau) -- you are of course quite correct -- it takes the equation into troubled waters. I skipped over it no doubt because the form CBro was using did not include the delay (another papers did not) -- but it was still a nonlinear system for which I suspected no closed form solution exists. I have to admit that I did not read the material in depth, thinking to comment only on the Mathematica aspects.

And thanks, Pat -- I have not used delays in diffeqs before. I knew generally that the capability exists, but I will look into it.

I don't know if Carolyn is still following this discussion, but I hope she rejoins it.

Best, David

Dave,

If it is more complicated that in the paper, then it is really complicated. I would settle for a good numerical solution.

I am not sure you understood my point: I was questioning whether the thing was a differential equation at all (linear or not). Once the behavior of your system depends on the value of some quantity some finite time ago, you lose locality in time. It may not be a differential equation even if it looks like one. For example, dx/dt = x(x(t)) is not a differential equation.

That may seem like a pedantic formality, but the reason I point it out is that if the thing is really not a differential equation, then most likely the approaches built into Mathematica do not work. Of course, there is always hope for a numerical approach.

As for Pittsburgh, what can I say, know I know two David Keiths. I am sure Vancouver is no worse...

OL

Hello Carolyn and David,

This looks like a really tough problem. I just saw this today and took a look at the paper, so I have not tried to solve it yet.

I have to say that if Mathematica cannot solve a differential equation symbolically, then you would be hard pressed to find another program that could. Not trying to say that Mathematica is the best thing ever (although I do think it is quite impressive), it is just that not many programs are good at symbolic manipulations. I recall sitting through a presentation about a MAPLE tool-bench that couples to MATLAB. They were discussing symbolic solution of differential equations in the context of space probe design or the like. You may want to check that. Unfortunately I do not have references handy.

I think part of the reason your problem becomes so hard is that when you introduce the functions evaluated at a different time in the right-hand side (x2(t-tau)), then strictly speaking you probably no longer have a differential equation (not any equation with derivatives in it is a differential equation). I will try to take a look when I have a chance, but it looks like a symbolic solution would be very hard.

I agree with David on the subscripts. They are nice to display results in elegant forms. I have spent hours chasing illogical results which at the end turn out to be a mess-up with a subscripted variable.

A separate question for David, did you use to live in Pittsburgh? It just so happens that I know a David Keith...

Best,

OL.

And here is a numerical solution for your nonlinear system, with made up parameter values. When forming systemB, I used Thread so equations would be formatted individually, so they could be joined with the initial conditions. I also attach the notebook.

In[1]:= B = {{(p0 - q0)*(v0/(1 + \[Beta]0 ((x2[t])^2))) - d0, 0, 
   0}, {(1 - p0 + 
      q0)*(v0/(1 + \[Beta]0 ((x2[t])^2))), (p1 - 
       q1)*(v1/(1 + \[Beta]1 ((x2[t])^2))) - d1, 
   0}, {0, (1 - p1 + q1)*(v1/(1 + \[Beta]1 ((x2[t])^2))), -d2}}

Out[1]= {{-d0 + ((p0 - q0) v0)/(1 + \[Beta]0 x2[t]^2), 0, 
  0}, {((1 - p0 + q0) v0)/(
  1 + \[Beta]0 x2[t]^2), -d1 + ((p1 - q1) v1)/(1 + \[Beta]1 x2[t]^2), 
  0}, {0, ((1 - p1 + q1) v1)/(1 + \[Beta]1 x2[t]^2), -d2}}

In[2]:= initB = {x0[0] == x00, x1[0] == x10, x2[0] == x20};

In[12]:= X[t_] = {x0[t], x1[t], x2[t]};

In[13]:= systemB = Thread[X'[t] == B.X[t]];

In[14]:= equations = Join[systemB, initB];

In[7]:= par = {q0 -> .5, q1 -> .5, v0 -> 1, v1 -> 1, x00 -> 100, 
   x10 -> 0, x20 -> 0, p0 -> .5, p1 -> .4, p2 -> .5, d0 -> .2, 
   d1 -> .3, d2 -> .4, \[Beta]0 -> 1, \[Beta]1 -> 1};

In[8]:= equations /. par

Out[8]= {Derivative[1][x0][t] == -0.2 x0[t], 
 Derivative[1][x1][t] == (1. x0[t])/(1 + x2[t]^2) + 
   x1[t] (-0.3 - 0.1/(1 + x2[t]^2)), 
 Derivative[1][x2][t] == -0.4 x2[t] + (1.1 x1[t])/(1 + x2[t]^2), 
 x0[0] == 100, x1[0] == 0, x2[0] == 0}

In[15]:= {fx0, fx1, fx2} = 
  NDSolveValue[equations /. par, {x0, x1, x2}, {t, 0, 100}];

In[17]:= LogPlot[{fx0[t], fx1[t], fx2[t]}, {t, 0, 100}, 
 PlotLegends -> {"x0", "x1", "x2"}, PlotTheme -> "Scientific", 
 FrameLabel -> {"Time(s)", "Concentration"}]