Thank you one and all for your helpful responses. It looks exactly like what I was looking for.

Maybe you want:

sol = DSolve[{d'[t] == C1 Exp[-g t] - c d[t], d[0] == D0}, d[t], t]

(*{{d[t] -> (E^(-c t) (-C1 + c D0 + C1 E^((c - g) t) - D0 g))/(c - g)}}*)

DSolve[{b'[t] == C2/(d[t] /. sol[[1]])^n - c b[t], b[0] == B0}, b[t], t] // Simplify

(*{{b[t] -> (1/(
   c (1 + n)))(B0 c E^(-c t) (1 + n) - 
     C2 D0^-n E^(-c t) (1 - C1/(C1 - c D0 + D0 g))^
      n Hypergeometric2F1[n, (c (1 + n))/(c - g), (-g + c (2 + n))/(
       c - g), C1/(C1 + D0 (-c + g))] + 
     C2 ((E^(-c t) (C1 (-1 + E^((c - g) t)) + D0 (c - g)))/(
       c - g))^-n (1 - (C1 E^((c - g) t))/(C1 - c D0 + D0 g))^
      n Hypergeometric2F1[n, (c (1 + n))/(c - g), (-g + c (2 + n))/(
       c - g), (C1 E^((c - g) t))/(C1 + D0 (-c + g))])}}*)

Or more shorter:

  ClearAll["`*"]; Remove["`*"];

  d[t] = DSolveValue[{d'[t] == C1 Exp[-g t] - c d[t], d[0] == D0}, d[t], t]
  DSolve[{b'[t] == C2/d[t]^n - c b[t], b[0] == B0}, b[t], t] // Simplify

More more shorter:

  ClearAll["`*"]; Remove["`*"];

  DSolve[{d'[t] == C1 Exp[-g t] - c d[t], d[0] == D0, 
     b'[t] == C2/d[t]^n - c b[t], b[0] == B0}, {d[t], b[t]}, t] // Simplify

   (*{{d[t] -> (E^(-c t) (C1 (-1 + E^((c - g) t)) + D0 (c - g)))/(c - g), 
    b[t] -> (1/(c (-C1 + D0 (c - g)) (1 + n)))
     E^(-c t) (B0 c (-C1 + D0 (c - g)) (1 + n) - 
        C2 D0^(1 - n) (c - g) Hypergeometric2F1[1, (2 c + g (-1 + n))/(
          c - g), 1 + (c (1 + n))/(c - g), C1/(C1 + D0 (-c + g))] + 
        C2 E^((c - g) t) (C1 (E^(c t) - E^(g t)) + 
           D0 E^(g t) (c - g)) ((
          E^(-c t) (C1 (-1 + E^((c - g) t)) + D0 (c - g)))/(
          c - g))^-n Hypergeometric2F1[1, (2 c + g (-1 + n))/(c - g), 
          1 + (c (1 + n))/(c - g), (C1 E^((c - g) t))/(
          C1 + D0 (-c + g))])}}*)

Thank you Daniel,

The equations are indeed coupled, which poses a new hardship in my attempts - putting aside the analytical solution (in which I used the parameters e [aka C1],g,n as fitted by another software plus initial conditions c, D0,B0), I set off to make the other software redundant by achieving the parameter fitting from Mathematica.

What I am able to achieve by now is fitting e,g over the d(t) ODE, then taking the fitted e,g and using them to fit n in the b(t) ODE.

However, my attempts at fitting e,g,n at a single stroke using the coupled equations for d(t) and b(t) have not been successful thus far. I am either not getting the syntax right, or possibly other commands should be used.

I wonder if this can be achieved in Mathematica (I would safely bet that the answer is yes)?

Here is how I fit in the "two-stroke" way (p1hdv and p1hbv are the datasets I am fitting against):

hdvodep = d'[t] == (1 - e1) c D0 Exp[-g1 t] - c d[t];

hdvmod = ParametricNDSolveValue[{hdvodep, d[0] == D0}, d, {t, 0, 28}, {e1, g1}];

fitd = FindFit[p1hdv, {Log[10, hdvmod[e1, g1][t]], e1 < 1.00, g1 > 0.01}, {e1, g1}, t]

which gives nice values for e1, g1 and then:

hdvsolp = DSolve[{hdvodep /. fitd, d[0] == D0}, d[t], t] // Simplify;

hbvodep = b'[t] == c B0 n1^-4 D0^n1 (d[t] /. hdvsolp[[1]])^-n1 - c b[t];

hbvmod = ParametricNDSolveValue[{hbvodep, b[0] == B0}, b, {t, 0, 28}, {n1}];

fitb = FindFit[p1hbv, {Log[10, hbvmod[n1][t]]}, {n1}, t]

which yields a fit for n1.