hello
try this
now beta and zeta must be vectors so that R_0 is a matrix (2D array)
use contourf for the plot
I let you choose which colormap fits best your color taste !
hope it helps
%Basic Reproduction Number
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
zeta = linspace(0,1,100);
% beta = 0.88; %transmission coefficient
beta = linspace(-0.001,1,100)'; % NB the min value must be slightly negative to let the contour line "0" show up
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
K_1 = eta_1 + eta_2 + eta_3 + mu;
K_2 = zeta + q + mu;
K_3 = gamma_1 + delta_1 + mu;
K_4 = alpha + gamma_2 + delta_2 + mu;
R_0 = beta.*(tau.*eta_2*K_2.*K_4 + K_3.*(eta_3.*K_2 + eta_1.*q))./(K_1.*K_2.*K_3.*K_4);
% plot
levels = (0:0.2:1.2);
[c,h] = contourf(beta,zeta,R_0',levels);
clabel(c,h)
h.LineWidth = 2; % contour line width (up to you)
xlabel('\beta');
ylabel('\zeta');
colormap(jet(numel(levels)-1));
colorbar('vert');
