How to verify if the stationary states of a biological model are correct? I'm mathematically modeling yellow fever, and yet, after assembling the model, I need to find the stationary states of the model. I had to do it by hand because wolfram software could not get the necessary values. 

My question is, how do I check if my stationary states are in line with my biological model? 

I was told that there is a way to verify this by performing a substitution of the s values of the points within my own template code. But I have no idea how to do it.

The code used is taking forever to get back to me as this data behaves in my model
    

function ypsir = ypsirtry(t,y)
betah = .01;
muh = 0.01;
sigmah = 0.05;
alphah = 0.06;
deltah = 0.05; 
mui = 0.07;
mum = 0.07;
betai = 0.03;
lambdai = 0.06;
lambdah = 0.5;
k = 500;

ypsir = zeros(5,1);
ypsir(1) = sigmah*y(3) - muh*y(1) - betah*y(1)*y(5) + lambdah*(y(2)+ y(3) + 
y(1));
ypsir(2) = -deltah*y(2) - y(2)*(alphah + muh) + betah*y(5)*y(1);
ypsir(3) = deltah*y(2) - sigmah*y(3);
ypsir(4) = -mui*y(4) - mum*y(4) - betai*y(2)*y(4) + lambdai*(y(4) + y(5))*(1 
- (y(4) + y(5))/k);
ypsir(5) = betai*y(2)*y(4) - mui*y(5) - mum*y(5);
end


Along with this calling program
    

clear;
to = 0;
tf =50;
yo = [99 1 0 1 3];
[t y] = ode45(@(t,y)ypsirtry(t,y),[to tf],yo); 
plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5))
title('First Model of Yellow Fever')
xlabel('time')
ylabel('susceptible, infected, recovered,Mosquito Non-carrier, Carrier 
mosquito') 


My stationary states are:
    

S = (1/betai)*[(muh+alphah)+(sigmah/deltah)*deltah]*(mui+mum)/betai*N*I);

I = sigmah*R/deltah;

R = ((muh/betah + lambdah)*(muh+alphah))/(sigmah+lamdah)-(muh+alphah)* 
(deltah/sigmah*deltah)-(betah/muh - lambdah)*(betai*N*I)* 
(delta/sigmah*deltah));

N = -k/lambdai((-lamdai+(lambdai/k))*P/-betai*I + lambdai - (2k/lambdai*P)- 
mui-mum);

P = betai*N*I/(mui+mum); 


At where
S = susceptible
I = infected
R = retweeted
N = non-carrier mosquito
P = carrier mosquito
After much analysis of the model, I came to the conclusion that I will not be able to verify how the equilibrium points behave graphically. Since you will have to create a phase diagram in R ^ 5, there are five equations, and I can not do two-to-two comparison in these cases, because I will always have a function that depends on a variable that will not be placed, but will be in the equation plotted anyway.

But in any case, if it were in three dimensions, for the purpose of information for those who are interested, the verification of the points would be to place the points (which in calculating have both numbers and not functions in function of other functions), within the model created, hence plot the phase diagram and check whether the arrows converge or not to the point. If they are converging, they are correct. 

Regarding the model is taking forever, I took the logistics part and ran well!