@dharr 

Thnak you for the idea. 
Using Routh-Hurwitz 'array'  can we solve the system to get the set of values where the roots are in the left half plane.

@mmcdara 

Thanks for your remarks. 

All parameters are positive in my equation. But, Maybe we can add  a condition so that we can determine or know the sign of the roots of the equation. 

My first goal is to study the stability of a given equilibrium point, so I computed the characteristic equation which contains  many paramters. So, if I know the sign of my roots I can know the nature of my equilibrium point. 

Maybe we can add some condition that gives the sign of the roots 

@mmcdara 

alpha>0, 
beta>0 
and 1>alpha+beta>0

and we suppose mu>1 

I tried also asymp  to get the asymptotic of F when t goes to infinity but unfortunately no result obtained

@acer 

I appreciate your comments

@Rouben Rostamian  

Thank you for your code. I really appreciate the result to verify my computation by hand. 
By, when I run the code, I obtain the following erros, I don't understand the kind of these errors.

Thank you

@tomleslie 

curve_three_d.mw

@Carl Love 

Thank you for your help

The modification proposed are added in the code,

But its seems that the x axis is used in the reverse order,

The  x-axis from 1740 to 4033 , but the figure show that the x-axis is  reversed,

@Carl Love 

Thank you

@tomleslie 

Thank you for your answer

@tomleslie 

My problem is solved based on your code 

thank you

@acer 

Thank you

 I appreciate your remarks

@tomleslie 

@tomleslie
According to all remarks received, I tried to do my best to correct the code.  

code1.mw

But wht the code does not add blocks of matrix K, 

This code, is exactly similar to matlab code presented in the following link, used to solve odes by finite element

https://www.youtube.com/watch?v=l4Y3iQrPBa8

thank you for all remarks 

@acer 

thank you for your remark

@acer 

compare other method with fsolve