# How to analytically solve for the impulse response function of an ODE system?

I'm quite new to Mathematica, and I'm quite rusty on my differential equations, so I apologise in advance if this is trivial!

I'm trying to use NDSolve to analytically solve a set of differential equations for a compartmental model over time (t) with an input function. Effectively, I want to figure out the impulse-response function of the following model:

 Input  |  System:  (CT = C1 + C2)
|
____  |  K1   ____    k3   ____
| C0 | |----> | C1 | ----> | C2 |
|____| |<---- |____| <---- |____|
|  k2           k4
|


where C0(t) is our known input function, measured over time (however for this exercise, I'm interested in C0(t) as if it were a Dirac delta function). We also measure the total of C1(t)+C2(t) = CT(t) (T for total). Hence, the convolution of our input function C0(t) and the impulse response function, H, generates our measured curve CT(t):

CT(t) = (H * Cp)(t)


The differential equations are as follows:

  dC1(t)/dt = K1*C0(t) - (k2+k3)*C1(t)
dC2(t)/dt = k3*C1(t) - k4*C2(t)


I'm not quite sure how to go about solving this in Mathematica. Ideally, I would like to simplify it to a sum of two exponential functions, as was done here (page 648, screenshotted below): Is there any way I do this? I've been trying the following:

NDSolve[{C1'[t] == K1*C0[t] - (k2 + k3)*C1[t] + k4*C2[t],
C2'[t] == k3*C1[t] - k4*C2[t],
C1 == 0,
C2 == 0}, {C1, C2}, {t, 0, 120}]


but I get the error:

NDSolve::underdet: There are more dependent variables, {C0[t],C1[t],C2[t]}, than equations, so the system is underdetermined.

I just want to be able to tell Mathematica that C0 is known, and that I only want C1 and C2, which can be convolved with the input function.

Any help would be greatly appreciated!

• What exactly is C0[t]? BTW for numerical solution you need to provide values for all k's.
– zhk
Feb 8 at 1:24
• You have 3 dependent variables: C0,C1,C2, but only 2 equations Feb 8 at 8:42
• Sorry for being unclear! C0[t] is the input to C1 and C2. In the experimental setting, I measure it over time and I fit a model to it, which I use as input. So it's a known and I'm not trying to solve for it. I'm trying to solve for H(t), the impulse response function, i.e. if C0(t) were a Dirac delta function, what would the function look like. For fitting the curve, I use nonlinear optimisation: that part is easy. I want to understand how I could use Mathematica to solve for Phi1, Phi2, Theta1 and Theta2 in the attached solution from the article. Feb 8 at 9:05
• NDSolve can't work with symbolic expressions (or symbolic functions). Maybe this is an indication of a direction Block[{C0, K1 = 1, k2 = 2, k3 = 3, k4 = 4}, C0[t_] := Exp[-t]; NDSolve[{C1'[t] == K1*C0[t] - (k2 + k3)*C1[t] + k4*C2[t], C2'[t] == k3*C1[t] - k4*C2[t], C1 == 0, C2 == 0}, {C1, C2}, {t, 0, 120}] ]  Feb 8 at 20:40
• Riight! I'm barking up the wrong tree with NDSolve then... I don't actually want to integrate through the function. I want it to solve for the symbolic impulse response function using Laplace transform. I think my question is ill-posed. Thank you so much: that's really helpful!! Feb 8 at 22:15

DSolve[{C1'[t] == K1*C0[t] - (k2 + k3)*C1[t] + k4*C2[t],    C2'[t] == k3*C1[t] - k4*C2[t], C1 == 0, C2 == 0}, {C1[t],    C2[t]}, t]