function ResponseToP
%% define parameters
p = 5;    % production rate of MnP per unit bacteria per day = 5 unit/min
b0 = 100; % starting concentration of bacteria               = 100 unit
v = 0.879;% volume                                           = 0.879 m^3
j = 0.19; % flow rate                                        = 0.18m^3/min

k1f = 2;  % forward reaction of Enzyme-substrate formation   = unit^-1 * min^-1
k1r = 0.5;% reverse reaction of Enzyme-substrate formation   = min^-1
k2 = 1;   % forward reaction of breaking down of nylon       = min^-1

n = 1000;   % Nylon concentration of waster water flowing into chamber = unit

params = {p,v,b0,j,k1f,k1r,k2};
%% Run simulation for a range of production rate, p
paramRange = [0.05;0.2;1;5;20];
% paramRange =[0;0.1;1;5;50];
MnPfinal = zeros(size(paramRange));
figure(3);
for i=1
    p = paramRange(i);
    params = {p,v,b0,j,k1f,k1r,k2};
    c0 = [0; p*b0; n; 0];
    tspan = [0 60];
    options = [];
    [t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
    subplot(1,2,1); plot(t,c(:,1),'b.-'); hold on;
    MnPfinal(i) = c(end,1);

end

for i=2
    p = paramRange(i);
    params = {p,v,b0,j,k1f,k1r,k2};
    c0 = [0; p*b0; n; 0];
    tspan = [0 60];
    options = [];
    [t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
    subplot(1,2,1); plot(t,c(:,1),'g.-'); hold on;
    MnPfinal(i) = c(end,1);

end

for i=3
    p = paramRange(i);
    params = {p,v,b0,j,k1f,k1r,k2};
    c0 = [0; p*b0; n; 0];
    tspan = [0 60];
    options = [];
    [t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
    subplot(1,2,1); plot(t,c(:,1),'r.-'); hold on;
    MnPfinal(i) = c(end,1);

end

for i=4
    p = paramRange(i);
    params = {p,v,b0,j,k1f,k1r,k2};
    c0 = [0; p*b0; n; 0];
    tspan = [0 60];
    options = [];
    [t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
    subplot(1,2,1); plot(t,c(:,1),'y.-'); hold on;
    MnPfinal(i) = c(end,1);

end
for i=5
    p = paramRange(i);
    params = {p,v,b0,j,k1f,k1r,k2};
    c0 = [0; p*b0; n; 0];
    tspan = [0 60];
    options = [];
    [t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
    subplot(1,2,1); plot(t,c(:,1),'c.-'); hold on;
    MnPfinal(i) = c(end,1);

end

% for i=6
%     p = paramRange(i);
%     params = {p,b0,v,j,k1f,k1r,k2};
%     c0 = [0; p*b0; n; 0];
%     tspan = [0 60];
%     options = [];
%     [t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
%     subplot(1,2,1); plot(t,c(:,1),'c.-'); hold on;
%     MnPfinal(i) = c(end,1);
% 
% end


hold off; 
% axis([-1 30 -0.05 1.05]);
legend('0.05','0.2','1','5','20')
xlabel('Time (hr)'); ylabel('Concentration of MnP (unit)');


paramRange = [0.05;0.1;0.2;0.5;1;5;20];
% paramRange =[0;0.1;1;5;50];
MnPfinal = zeros(size(paramRange));

for i=1:length(paramRange)
    p = paramRange(i);
    params = {p,v,b0,j,k1f,k1r,k2};
    c0 = [0; p*b0; n; 0];
    tspan = [0 60];
    options = [];
    [t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
  
    MnPfinal(i) = c(end,1);
end
subplot(1,2,2)
plot(paramRange,MnPfinal,'bo-');
xlabel('Production rate of MnP per unit number bacteria (unit/number of bacteria)'); ylabel('Concentration of MnP after 60 hours(unit)');
% legend('MnP');


