Setting assumptions for a function when solving system of ODEs

I have a following code below that is supposed to solve a dynamic system. I am trying to restrict X3(t) to be less than or equal to 0.1 (a saturation point). So I added a line about assumption for X3(t). However, it seems like assume command is not possible for symbolic functions. How should I get this done?

syms t A B kEuT kTEu X1(t) X2(t) X3(t) X4(t) X5(t) X6(t)

eqn1 = diff(X1) == -B*X1 + A*X2
eqn2 = diff(X2) == -(A+kEuT)*X2 + B*X1 + kTEu*X3
eqn3 = diff(X3) == kEuT*X2 - kTEu*X3
C1 = X1(-20) == 1
C2 = X2(-20) == 0
C3 = X3(-20) == 0
assume(X3 <= 0.1)
[X1S(t) X2S(t) X3S(t)] = dsolve(eqn1, eqn2, eqn3, C1, C2, C3)
X1Sc = simplify(combine(X1S(t)),'IgnoreAnalyticConstraints',true)
X2Sc = simplify(combine(X2S(t)),'IgnoreAnalyticConstraints',true)
X3Sc = simplify(combine(X3S(t)),'IgnoreAnalyticConstraints',true)

eqn4 = diff(X4) == A*X5
eqn5 = diff(X5) == -(A+kEuT)*X5 + kTEu*X6
eqn6 = diff(X6) == kEuT*X5 - kTEu*X6
C4 = X4(0) == X1S(0)
C5 = X5(0) == X2S(0)
C6 = X6(0) == X3S(0)
[X4S(t) X5S(t) X6S(t)] = dsolve(eqn4, eqn5, eqn6, C4, C5, C6)
X4Sc = simplify(combine(X4S(t)),'IgnoreAnalyticConstraints',true)
X5Sc = simplify(combine(X5S(t)),'IgnoreAnalyticConstraints',true)
X6Sc = simplify(combine(X6S(t)),'IgnoreAnalyticConstraints',true)

X1Sc = vpa(subs(X1Sc, [A, B, kEuT, kTEu], [2000000, 1000000, 20000, 100]))
X2Sc = vpa(subs(X2Sc, [A, B, kEuT, kTEu], [2000000, 1000000, 20000, 100]))
X3Sc = vpa(subs(X3Sc, [A, B, kEuT, kTEu], [2000000, 1000000, 20000, 100]))
X4Sc = vpa(subs(X4Sc, [A, B, kEuT, kTEu], [2000000, 1000000, 20000, 100]))
X5Sc = vpa(subs(X5Sc, [A, B, kEuT, kTEu], [2000000, 1000000, 20000, 100]))
X6Sc = vpa(subs(X6Sc, [A, B, kEuT, kTEu], [2000000, 1000000, 20000, 100]))

x1 = linspace(-20,0,200)
X1Scp = subs(X1Sc, t, x1)
X2Scp = subs(X2Sc, t, x1)
X3Scp = subs(X3Sc, t, x1)
x2 = linspace(0,20,200)
X4Scp = subs(X4Sc, t, x2)
X5Scp = subs(X5Sc, t, x2)
X6Scp = subs(X6Sc, t, x2)

plot(x1,X1Scp,'b',x1,X2Scp,'k',x1,X3Scp','r',x2,X4Scp,'b',x2,X5Scp,'k',x2,X6Scp,'r'), xlabel('time/s'),ylabel('population density'),axis([-20 20 0 0.01])

Thank you

Answer 1

[Nicolas Mira Gebauer]

Hi, how did you solve this problem?

I am currently stuck at something similar..

Thank you in advance