@Preben Alsholm 

It makes some sense. Thank you so much for your wonderful contribution. I will surely require your help again soon, I am still getting myself equiped with maple and you have made it seem so easy to grasp with your help. Thank you!!!

@Preben Alsholm 

Thank you so much for the help. I am still getting my head around the elimination part where you suggest to get rid off the entire term....

(b*(I1(t)+phi1*I2(t)+phi2*I3(t))/N(t))...when N(t)->infinity as t->infinity.

The thing is although it makes a difference in doing so, I am still wondering if there isnt an alternative way without getting rid off this term completely....as I think it has significant influence especially in the solution set or equilibria? Just wondering!

@Preben Alsholm 

Yes it is irrelevant, but is there a way to get rid off the negative values and obtain only positive values?

Thank you in anticipation.

@Preben Alsholm 

 Yes it worked but i still get a weird equilibria coordinates (strange solution) when i use the following,

##Equilibria

solve{rhs~({sys1,sys2,sys3,sys4,sys5,sys5,sys6}),{S(t),I1(t), I2(t), I3(t) ,I4(t), A(t)});

eval([%],params);

Map(subs,%[ S(t),I1(t), I2(t), I3(t) ,I4(t), A(t)}]);

And I dont seem to see any mistake. please help.

Thanks in ancticipation.

@THAPELO

Hi there,

Thank you so much I tried your methods on three and four equations and it works howerver I don't know where I am wrong when I try it on six equations.....

This is the ODE system i am working on:

@Preben Alsholm 

Thank you so much for sharing the light :). I will definitely try this method on my six ode system.

@Preben Alsholm

Thank you so much for your comment. Please excuse me as I am still new in maple community(typing* etc). This is what I obtain, when I put numeric at the end as you suggests,

.....and still no results!!

This is the system in detail

The original system is as follows

dS1/dt = b - d*S1 - (k*S1*S2)/(1+ alpha*S1) - gamma*S3,

dS2/dt = (k*S1*S2)/(1+ alpha*S1) - (d + mu)*S2,

dS3/dt = mu*S2 - (d + gamma)*S3.

The parameters are taken as b=1, d=0.2, k=1, alpha=0.5, gamma=0.25 and mu=0.3.

The initial conditions are taken as S1(0)=1.35, S2(0)=0.9, S3(0)=0.45.

sys1:= (diff(S1(t),t) = b - d*S1(t) - (k*S1(t)*S2(t))/(1+ alpha*S1(t)) - gamma*S3(t))

sys2:= (diff(S2(t),t) = (k*S1(t)*S2(t))/(1+ alpha*S1(t)) - (d + mu)*S2(t))

sys3:= (diff(S3(t),t) = mu*S2(t) - (d + gamma)*S3(t))

Sol1:=dsolve({sys1, sys2, sys3, S1(0)=1.35, S2(0)=0.9, S3(0)=0.45}, [S1(t), S2(t), S3(t)], numeric)

This is what i get: Error, (in dsolve)  too many arguments; some or all of the following are wrong: [{S1(t), S2(t), S3(t)}, [S1(t), S2*t, S3(t)]]

All i need is to obtain numeric values of S1(t), S2(t) and S3(t).

Thank you in anticipation.