Dear all, this post arises from my other post in Chemical Engineering about integration of kinetic equations.

I was wondering if there is any way to force or restrict integration via NDSolve to only give answers in the real domain. At the minute, I am producing complex answers to real equations which is rather tricky to work with! I have seen a few posts around this topic but haven't come across a comprehensive answer yet. I was wondering if anyone with more experience than me may be able to offer their help.

Many thanks! Robin

P.S. I must say a big thank you to David Keith who has been incredibly helpful with this matter so far!

The code I am running is listed below.

Part 1:

ClearAll[CSA, CAA, CASA, CHA, CASAA, CH2O, V, t]

equs = {CAA'[t] == -(k1*CSA[t]*CAA[t]) - (k2*CASA[t]*CAA[t]) - (1/2)*(k4*CAA[t]*CH2O[t]) - CAA[t]*((V'[t])/(V[t])),
CHA'[t] == (k1*CSA[t]*CAA[t]) + (k2*CASA[t]*CAA[t]) + (k3*CASAA[t]*CH2O[t]) + (k4*CAA[t]*CH2O[t]) - CHA[t]*((V'[t])/(V[t])),
CSA'[t] == (kd*(CSatSA - CSA[t])^d) - (k1*CSA[t]*CAA[t]) - CSA[t]*((V'[t])/(V[t])),
CASAA'[t] == (k2*CASA[t]*CAA[t]) - (k3*CASAA[t]*CH2O[t]) - CASAA[t]*((V'[t])/(V[t])),
CH2O'[t] == -(k3*CASAA[t]*CH2O[t]) - (1/2)*(k4*CAA[t]*CH2O[t]) - CH2O[t]*((V'[t])/(V[t])),
CASA'[t] == (k1*CSA[t]*CAA[t]) - (k2*CASA[t]*CAA[t]) + (k3*CASAA[t]*CH2O[t]) - (kc*(CASA[t] - CSatASA)^c) - CASA[t]*((V'[t])/(V[t])),
V'[t] == (V[t]*(vASA*((k1*CSA[t]*CAA[t]) - (k2*CASA[t]*CAA[t]) + (k3*CASAA[t]*CH2O[t]) - (kc*(CASA[t] - CSatASA)^c)) + 
   vH2O*(-(k3*CASAA[t]*CH2O[t]) - (1/2)*(k4*CAA[t]*CH2O[t]) + CInH2O*(f/V[t])) + vASAA*((k2*CASA[t]*CAA[t]) - (k3*CASAA[t]*CH2O[t])) + 
   vSA*((kd*(CSatSA - CSA[t])^d) - (k1*CSA[t]*CAA[t])) + vHA*((k1*CSA[t]*CAA[t]) + (k2*CASA[t]*CAA[t]) + (k3*CASAA[t]*CH2O[t]) +  
  (k4*CAA[t]*CH2O[t])) + vAA*(-(k1*CSA[t]*CAA[t]) - (k2*CASA[t]*CAA[t]) - (1/2)*(k4*CAA[t]*CH2O[t]))))}

initial = {CAA[0] == (1057/100), CHA[0] == 0, CSA[0] == 0, CASAA[0] == 0, CH2O[0] == 0, CASA[0] == 0, V[0] == 120}

vals = {k1 -> (337/10000), k2 -> (43/100), k3 -> 1525, k4 -> 85, kd -> (44/100), kc -> (109/1000), c -> (134/100), d -> (190/100), 
  CSatASA -> (152/100), CSatSA -> (219/100), vASA -> (14/1000), vSA -> (14/1000), vASAA -> (14/1000), vHA -> (14/1000), vAA -> (14/1000), 
  vH2O -> (14/1000), CInH20 -> (555/10), f -> 0}

Part 2:

{sCSA, sCASA, sCASAA, sCHA, sCH2O, sCAA, sV} = Reaction = NDSolveValue[{initial, equs /. vals}, {CSA[t], CASA[t], CASAA[t], 
 CHA[t], CH2O[t], CAA[t], V[t]}, {t, 0, 200}];

myplot = Plot[Reaction // Evaluate, {t, 0, 200}, PlotPoints -> 500, PlotRange -> {{0, 200}, {0, 20}}, AxesLabel -> {time, conc}, 
LabelStyle -> Directive[Black, Medium, "Palatino"], 
  PlotLegends -> {CSA, CASA, CASAA, CHA, CH2O, CAA, V}]