*I have been working on drawing the contour plots of the basic reproduction number (R0) of a model with respect to two parameters. For practice, I am currently working on replicating the graph attached below from a paper by Kifle et al,

The MATLAB codes that I used to calculate R0 of the model are given below.

%Parameters 

theta = 141302; %recruitment rate
mu = 0.001229; %natural death rate
tau = 0.45; %modification factor for A
zeta = 1/14; %influx from Q to S 
beta = 0.88; %transmission coefficient
alpha = 0.75214; %hospitalization rate
q = 0.31167; %influx from Q to I
eta_1 = 0.81692; %influx from E to Q
eta_2 = 0.02557; %influx from E to A
eta_3 = 1/7; %influx from E to I
delta_1 = 0.16673; %disease death rate for A
delta_2 = 0.00147; %disease death rate for I
delta_3 = 0.00038; %disease death rate for J
gamma_1 = 0.00827; %recovery rate for A
gamma_2 = 0.00787; %recovery rate for I
gamma_3 = 0.20186; %recovery rate for J