# Slow Solution from DSolve when Adding Initial Conditions to system of ODEs

I am trying to find an analytical solution to the following system of ODEs:

uB'[t] == -a*uB[t]*uC[t] + b*uD[t],
uC'[t] == -a*uB[t]*uC[t] + b*uD[t],
uD'[t] == a*uB[t]*uC[t]-b*uD[t]


I can find the numerical solution readily with the following:

Assuming[uB > 0 && uC > 0 && uD > 0,
NumericalSol = NDSolve[{
uB'[t] == -uB[t] uC[t] + uD[t],
uC'[t] == -uB[t] uC[t] + uD[t],
uD'[t] == uB[t] uC[t] - uD[t] ,
uB[0] == 1,
uC[0] == 1,
uD[0] == 0},
{uB[t], uC[t], uD[t]}, {t, 0, 1}]];


And I can find the general closed-form solution with this:

ExactSol = DSolve[{
uB'[t] == -a*uB[t] uC[t] + b*uD[t],
uC'[t] == -a*uB[t] uC[t] + b*uD[t],
uD'[t] == a*uB[t] uC[t] - b*uD[t]},
{uB[t], uC[t], uD[t]}, t] // Simplify


But as soon as I include initial conditions to my analytical solution to reproduce the numerical solution:

ExactSol = DSolve[{
uB'[t] == -a*uB[t] uC[t] + b*uD[t],
uC'[t] == -a*uB[t] uC[t] + b*uD[t],
uD'[t] == a*uB[t] uC[t] - b*uD[t],
uB[0] == 1,
uC[0] == 1,
uD[0] == 0},
{uB[t], uC[t], uD[t]}, t] // Simplify


The code tries to evaluate but never finishes - is there a clear reason why the constrained solution would be difficult to obtain? What am I missing?

Thanks for any help.

it is having hard time solving for $$c_i$$, from initial conditions. We can do it manually

ode1 = uB'[t] == -a*uB[t] *uC[t] + b*uD[t];
ode2 = uC'[t] == -a*uB[t] *uC[t] + b*uD[t];
ode3 = uD'[t] == a*uB[t] *uC[t] - b*uD[t];
ic = {uB[0] == 1, uC[0] == 1, uD[0] == 0};
exactSol =
First@DSolve[{ode1, ode2, ode3}, {uB[t], uC[t], uD[t]}, t] //
FullSimplify


eq1 = uB[t] /. exactSol
eq2 = uC[t] /. exactSol
eq3 = uD[t] /. exactSol


sol = Solve[{1 == eq1 /. t -> 0, 1 == eq2 /. t -> 0,
0 == eq3 /. t -> 0}, {C[1], C[2], C[3]}]


 sol = First@Normal[sol, ConditionalExpression]


 exactSol /. sol


Notice, the solutions are not unique!. There is arbitrary $$c_4$$ in them. You set this to any value in Z (I take this back. Any integer value for $$c_4$$ gives the same solution. So any integer value for C[4] will work).

 exactSol /. sol/. C[4]->1


## Verify against numerical

Assuming[uB > 0 && uC > 0 && uD > 0,
NumericalSol =
First@NDSolve[{uB'[t] == -uB[t] uC[t] + uD[t],
uC'[t] == -uB[t] uC[t] + uD[t], uD'[t] == uB[t] uC[t] - uD[t],
uB[0] == 1, uC[0] == 1, uD[0] == 0}, {uB, uC, uD}, {t, 0, 1}]];
Plot[Evaluate[{uB[t], uC[t], uD[t]} /. NumericalSol], {t, 0, 1}]


newSol = exactSol /. sol /. C[4] -> 1;
Plot[Evaluate[{uB[t], uC[t], uD[t]} /. newSol /. {a -> 1,
b -> 1}], {t, 0, 1}]