Drawing a heatmap of the max value of a state variable against two parameters in the system

I am working with an epidemic ODE model and I wish to draw a heatmap of y(4) Max against the two parameters p and eta. Any help would be greatly appreciated. Thank you!!!

function dy=EpiModel(t,y)  

dy = zeros(5,1); 

%Model paramters 

c = 0.001;
theta = 0.1134;
beta = 0.4;
p = 0.05;
eta = 0.05;
epsilon = 0.084;
gamma = 0.2;
delta = 0.011;
mu = 0.000011;

%Equations

dy(1) = theta - (p+c*beta*y(4)+mu)*y(1);

dy(2) = p*y(1) - (eta*c*beta*y(4)+mu)*y(2);

dy(3) = c*beta*y(4)*y(1) + eta*c*beta*y(4)*y(2) - (epsilon+mu)*y(3);

dy(4) = epsilon*y(3) - (gamma+delta+mu)*y(4);

dy(5) = gamma*y(4) - mu*y(5);

end

[Bas123]

Does someone have any inputs? Thank you!

[Torsten]

For each pair (p,eta), you get a function for y(4) over the time interval you integrate. So what do you want to plot ? Is there a stationary value for y(4) it converges to which should be the value for y(4) in the heat map ?

[Bas123]

Hi, thanks for your question. It is the maximum value of y(4), denoted, (i.e., peak of infections) that I want to plot with respect to the variation of p and η.

[Torsten]

And what are your initial values for y1-y5 ?

[Bas123]

Oh sorry, I should have specified them earlier. Here are the initial values: .

Answer 1

[Torsten]

tspan = 0:0.01:100;
y0 = [100000;0;0;1;0];
p = 0.001:0.001:0.01;
eta = 0.001:0.001:0.02;
for i = 1:numel(p)
  for j = 1:numel(eta)
    [~,Y] = ode45(@(t,y)EpiModel(t,y,p(i),eta(j)),tspan,y0);  
    maxy4(i,j) = max(Y(:,4));
  end
end
contourf(eta,p,maxy4)
colorbar

function dy=EpiModel(t,y,p,eta)  
...
end

[Bas123]

This is great, thank you so much!