@Carl Love 

Hi Carl, I've tried and succeeded in changing the symbol for imaginary unit. Great.

@Preben Alsholm 

ok. I didn't notice that I have 'I' included in the name of one of my unkowns, which may be considered as the imaginary part.

Thanks, problem solved. It would be great if Maple could allow us to declare 'I' not to be the imagainary part.

@Preben Alsholm 

Thanks Preben, it's great.     By doing eval(sys2,S[b, 1, O](t) = 2), the system could be simplified significantly, from a DAE to ODE. :))) This will definitely work. And , I really like the way you obtain all the unknowns from odes (convert(indets(rhs~(odes),function(identical(t))),list);), I've never tried this way before, but I think I will pick up this fuction and apply it in the future.

and thanks again for updating your answer by adding addtional approach on using piecewise functions.

Preben, could you please maybe try to keep the system as a DAE, of course all the obvious mistakes for instance division by zero must be eliminated first. I hope to keep the system as a DAE system, and solve it without doing eval(sys2,S[b, 1, O](t) = 2) . The reason why I insisted on doing so is that this is the simplified version of my model, and soon, I will expand my model and introduce serveral more algebric equations into the system. 

@Preben Alsholm 

btw, it maybe possible to convert equation S[b, 1, O](t) = 2 to a variable S[b,1,O]=2 and substitute this variable in all my differential equations, but I simply want to include this algebraic euqation in my DAE system, and solve it as a DAE SYSTEM.

@Preben Alsholm 

yes, here are the equations UNION initial values, there is only one algebraic equation (see below in bold and underlined).

{diff(S[b, 1, I](t), t) = 2960.040000-(74001/1000)*S[b, 1, I](t), diff(S[b, 1, ALK](t), t) = 222003/500-(74001/1000)*S[b, 1, ALK](t)-0.2285714286e-1*S[b, 1, S](t)*S[b, 1, O](t)*X[b, 1, B, H](t)/((10.0+S[b, 1, S](t))*(.2+S[b, 1, O](t)))+0.2764160564e-1*S[b, 1, S](t)*X[b, 1, B, H](t)/((10.0+S[b, 1, S](t))*(.2+S[b, 1, O](t)))-.3004761905*S[b, 1, NH](t)*S[b, 1, O](t)*X[b, 1, B, A](t)/((1+S[b, 1, NH](t))*(.4+S[b, 1, O](t)))+0.3571428571e-2*S[b, 1, ND](t)*X[b, 1, B, H](t), diff(S[b, 1, ND](t), t) = 747.4101000-(74001/1000)*S[b, 1, ND](t)-0.5e-1*S[b, 1, ND](t)*X[b, 1, B, H](t)+(3*(S[b, 1, O](t)/(.2+S[b, 1, O](t))+.3200000000/(.2+S[b, 1, O](t))))*X[b, 1, ND](t)/(.1+X[b, 1, S](t)/X[b, 1, B, H](t)), diff(S[b, 1, NH](t), t) = 925.0125000-(74001/1000)*S[b, 1, NH](t)-.320*S[b, 1, S](t)*S[b, 1, O](t)*X[b, 1, B, H](t)/((10.0+S[b, 1, S](t))*(.2+S[b, 1, O](t)))-.1024000000*S[b, 1, S](t)*X[b, 1, B, H](t)/((10.0+S[b, 1, S](t))*(.2+S[b, 1, O](t)))-2.123333334*S[b, 1, NH](t)*S[b, 1, O](t)*X[b, 1, B, A](t)/((1+S[b, 1, NH](t))*(.4+S[b, 1, O](t)))+0.5e-1*S[b, 1, ND](t)*X[b, 1, B, H](t), diff(S[b, 1, NO](t), t) = -(74001/1000)*S[b, 1, NO](t)-.2204362801*S[b, 1, S](t)*X[b, 1, B, H](t)/((10.0+S[b, 1, S](t))*(.2+S[b, 1, O](t)))+2.083333334*S[b, 1, NH](t)*S[b, 1, O](t)*X[b, 1, B, A](t)/((1+S[b, 1, NH](t))*(.4+S[b, 1, O](t))), diff(S[b, 1, S](t), t) = 592008/125-(74001/1000)*S[b, 1, S](t)-5.970149252*S[b, 1, S](t)*S[b, 1, O](t)*X[b, 1, B, H](t)/((10.0+S[b, 1, S](t))*(.2+S[b, 1, O](t)))+1.910447761*S[b, 1, S](t)*X[b, 1, B, H](t)/((10.0+S[b, 1, S](t))*(.2+S[b, 1, O](t)))+3*X[b, 1, S](t)*(S[b, 1, O](t)/(.2+S[b, 1, O](t))+.3200000000/(.2+S[b, 1, O](t)))/(.1+X[b, 1, S](t)/X[b, 1, B, H](t)), diff(X[b, 1, I](t), t) = 4.444444444-(1/9)*X[b, 1, I](t), diff(X[b, 1, ND](t), t) = 2.033333333-(1/9)*X[b, 1, ND](t)+0.2256e-1*X[b, 1, B, H](t)+0.3760e-2*X[b, 1, B, A](t)-(3*(S[b, 1, O](t)/(.2+S[b, 1, O](t))+.3200000000/(.2+S[b, 1, O](t))))*X[b, 1, ND](t)/(.1+X[b, 1, S](t)/X[b, 1, B, H](t)), diff(X[b, 1, P](t), t) = -(1/9)*X[b, 1, P](t)+0.24e-1*X[b, 1, B, H](t)+0.40e-2*X[b, 1, B, A](t), diff(X[b, 1, S](t), t) = 160/9-(1/9)*X[b, 1, S](t)+.276*X[b, 1, B, H](t)+0.460e-1*X[b, 1, B, A](t), diff(X[b, 1, B, A](t), t) = 0.1111111111e-3-.1611111111*X[b, 1, B, A](t)+.5*S[b, 1, NH](t)*S[b, 1, O](t)*X[b, 1, B, A](t)/((1+S[b, 1, NH](t))*(.4+S[b, 1, O](t))), diff(X[b, 1, B, H](t), t) = 32/3-.4111111111*X[b, 1, B, H](t)+4.0*S[b, 1, S](t)*S[b, 1, O](t)*X[b, 1, B, H](t)/((10.0+S[b, 1, S](t))*(.2+S[b, 1, O](t)))+1.280000000*S[b, 1, S](t)*X[b, 1, B, H](t)/((10.0+S[b, 1, S](t))*(.2+S[b, 1, O](t)))-3*X[b, 1, S](t)*(S[b, 1, O](t)/(.2+S[b, 1, O](t))+.3200000000/(.2+S[b, 1, O](t)))/(.1+X[b, 1, S](t)/X[b, 1, B, H](t)), S[b, 1, I](0) = 0, S[b, 1, ALK](0) = 0, S[b, 1, ND](0) = 0, S[b, 1, NH](0) = 0, S[b, 1, NO](0) = 0, S[b, 1, O](t) = 2, S[b, 1, S](0) = 0, X[b, 1, I](0) = 0, X[b, 1, ND](0) = 0, X[b, 1, P](0) = 0, X[b, 1, S](0) = 0, X[b, 1, B, A](0) = 0, X[b, 1, B, H](0) = 0}

@Markiyan Hirnyk 

Hi Markiya,

Thanks for your answer, and it is very kind of you to attach the code in the attachment.

I managed to solve my answer in a good way, and your answer is of great help.

@Preben Alsholm 

Thanks Preben,

the general idea of the solution to this problem is to use piecewise function, and with your answer, i may understand the correct way to call piecewise, and my problem is solved.

exactly, it is a little bit 'surprising' that Y(t) did not show a graph that a piecewise defined function should have, maybe it is, but not that obvious, the ' switch ' point is not so obvious on the figure.