% Created by Nick Fidalgo, Chris Gatesman, Phillip Less, Brett Peters, Samantha Sandwick, and Alex Shetler
% Hill function parameters
n = 50;
k_m = linspace(1*10^-5,.0015,n);
k_5 = linspace(30,65,n);
[K_5, K_m] = meshgrid(k_5, k_m);
om = 1;
% Parameters
sig = 5;
gam = 483;
phi = .35;
nu = 1;
mu = .2;
% Initial Conditions
B_0 = 1;
D_0 = 0;
E_0 = 1;
F_0 = 1;
xi = 1.8741*10^-4;

tspan = linspace(0,21);
S0 = [B_0;(D_0/B_0);E_0];
S_0 = [B_0;(D_0/B_0);E_0;F_0];
% A(t) parameters and dosages
%P = [.3 .6 1.2 2.4 3.6 4.8; 11.29 14.29 22.44 69.44 76.54 110.66; .82 .57 .58 .2 .3 .17];
%[m,n] = size(P);
A_0 = 76.54;
a = .3;

for i = 1:n
    for j = 1:n
        [t,H] = ode45(@(t,s) hill_M2(t,s,K_m(i,j),K_5(i,j),om,sig,gam,xi,phi,A_0,a),tspan,S0);
        E(i,j) = max(H(:,3));
        %E_n(i,j) = max(H_n(:,4));
    end
end
V = [1.5 ,2, 3, 4, 5];
%surf(K_5,K_m, E)
contour(K_5,K_m,E, V, 'ShowText', 'on','Linewidth',1.3)
xlabel('K_5_0')
ylabel('K_m_a_x')
title('Contour Map of E_m_a_x (3.6 mg/kg)')

%zlabel('E max')