@Rouben Rostamian  I am sorry I did not get back to you as regards your last suggestion. I had some personal issues and I had to stop my research for a while. However, I am back now. I  have tried carrying out a stability analysis using Linearization. However, I am still having some issues. My aim is to have a Theorem in my research stating that the Disease free equilibrium is locally asymptotically stable if Rvm (Basic reproduction number) < 1, and that the endemic equilibrium is locally asymptotically stable if Rvm (Basic reproduction number) > 1.

Below is my Jacobian matrix

                            ROW 1                           ROW 2                       ROW 3                         ROW 4

Column 1:=                                                                                                      

Column 2:=                 -psi                                     0                                   M

Column 3:=                                                                                                           

Column 4:=                 0                                0                                0                        +  M +  -

From the above matrix, I found out two things:

1. From the above Jacobian matrix, A negative eigenvalue is not obtainable,
2. By using the routhz hurwitz criterion, i reduced the matrix block from 4 by 4 to 3 by 3 by using column 4 as co-factor, however, after this matrix reduction,i am unable to get all terms relating to my basic reproduction number for me to compare and see if the conditions Rvm < 1 holds.

My Basic Reproduction Number is as follows:

Are there any suggestions you can give me to solve my problem?

@Rouben Rostamian  

Thanks. I will check it out.

@Rouben Rostamian  

To determine whether or not equilibrium points (Disease-Free or Endemic) of a nonlinear ordinary differential equation is globally stable, a Lyapunov function is often employed. Since there are no general methods for constructing Lyapunov functions, May I know if MAPLE can be used to determine these Lyapunov functions? I have been on it for a while but I do not seem to understand how other authors create their Lyapunov functions from the respective journals I have read. Thanks

@Rouben Rostamian  

Thank for your comments and your reply to my post.

For the first comment, it's not as if I was dismissing the solution when I=0. My advisor was much interested in the second solution. Computing by hand was very difficult that is why I posted the question initially. 

As regards the stability analysis, I am aware of the Routh-Hurwitz theorem and I will do as you have advised.

Happy New year to you and thanks once again for all your comments and suggestions.

@Rouben Rostamian  

 Thank you for taking out time to provide me with solutions to my question.

I appreciate all your efforts.

I am currently working on the file you uploaded and I think I should be able to make a conclusion as regards the sign of the coefficients.

In case I have any question, I will sure let you know.

Thanks so much for your time.

By the way, can disease-free equilibrium be obtained in a similar way as the endemic equilibrium point? Am just asking to satisfy my curiosity. I used Gaussian elimination method and will just like to know if there are other ways to obtain the disease-free for future purpose.

@Rouben Rostamian  

Thank you for taking out time to provide me with a solution to the question I posted.

As regards me checking the model, you are right. In my model, i have 5 compartments. However, since the last compartment isn't involved in all other compartments, I decided not to use in my analysis. When the last compartment is added together with the other compartments, we then obtain a constant.

I have carried out the steps as you have suggested and I obtained a solution. However, the solution is not in a format useful to me. I hoped to obtain a polynomial solution to the "I^th" term then I can use the Descartes rule of sign change to find out whether an endemic equilibrium point exists or not. 

Computing by hand to obtain the Ith term in a polynomial format is cumbersome and hectic. Is there any other way or method you can suggest to me?

Thanks for your help

@Yee Voon 

Can you please explain more? Where should the command be added? Thanks

@Rouben Rostamian  

I have tried using the idea you posted to solve my system of equations but I was unable to obtain a solution, My major concern is the "I " term. Since at endemic equilibrium state, "I is not equal to zero", how do i go about this? Tried solving by hand without any success.

Here is the system. Thanks

fS := (1-rho__M- rho__V)*Lambda-I*beta*S+omega__M*M+omega__V*V-(gamma__M+mu)*S;

fM := Lambda*rho__m+gamma__M*S-(I*beta*epsilon__M+mu+omega__M+varphi)*M;

fV := Lambda*rho__V+gamma__V*S+varphi*M-(I*beta*epsilon__V+mu+omega__V)*V;

fI := (M*beta*epsilon__M+V*beta*epsilon__V+S*beta-kappa-mu-tau)*I