@Carl Love @tomleslie

Thank you very much for your help, and I am very sorry for the late response, I am a new member of maple prime and my account was inactive for awhile.

The code works  and I manage to get some results. I also applied different models to my ODEs   and I end up having  expressions that  contain derivatives of the same order as the highest-order derivatives in the ODEs and I am working on it. 

@Carl Love, I am just wondering,  how can one get the  numeric integration of expressions of the dependent variables?  for example  if I need to integrate (\phi(t) (1+t) dt) and the limits of the integration for example are from 0 to t.  How is this possible? 

@Carl Love 

Thank you very much for your help, I have been stuck on this for almost a week. I have tried your code and I do not know why for the plotting part where
plots:-odeplot(dsol, eval[recurse](`~`[`[]`](t, ([seq])(`~`[rhs](Extra))), {Param_vals[], seq(Extra)}), t = 0 .. 2, legend = ([seq])(`~`[lhs](Extra)));

it gives me error
Error, invalid input: rhs expects 1 argument, but received 4. 

So instead of (Extra), I tried to use only one equation for example( Eq. 187), and see if I will get something different  but unfortunately, I still have the same problem