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):

Here is the solution I'd like to get with Mathematica

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] == 0,
         C2[0] == 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!

  • $\begingroup$ What exactly is C0[t]? BTW for numerical solution you need to provide values for all k's. $\endgroup$
    – zhk
    Feb 8, 2023 at 1:24
  • $\begingroup$ You have 3 dependent variables: C0,C1,C2, but only 2 equations $\endgroup$ Feb 8, 2023 at 8:42
  • $\begingroup$ 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. $\endgroup$
    – Granville
    Feb 8, 2023 at 9:05
  • $\begingroup$ 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] == 0, C2[0] == 0}, {C1, C2}, {t, 0, 120}] ] $\endgroup$ Feb 8, 2023 at 20:40
  • $\begingroup$ 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!! $\endgroup$
    – Granville
    Feb 8, 2023 at 22:15

1 Answer 1


DSolve will give you a symbolic result, albeit messy and containing integrals containing the arbitrary function C0[t]:

DSolve[{C1'[t] == K1*C0[t] - (k2 + k3)*C1[t] + k4*C2[t],    C2'[t] == k3*C1[t] - k4*C2[t], C1[0] == 0, C2[0] == 0}, {C1[t],    C2[t]}, t]
  • $\begingroup$ This is just what I was looking for! Thanks so much! $\endgroup$
    – Granville
    Feb 9, 2023 at 10:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.