% function returns best case equilibrium profit (vector over all T1 and T2 provided)
function [pmax,bw1,bc1] = max_profit2(t0,tend,m1,m2,T1,T2,R,bw2,bc2,min_bw,min_bc,Nd,Nw,Nc,IC,initial_guess)
% vector of best base pays
bw1 = zeros(length(T2),1); bc1 = zeros(length(T2),1); pmax = zeros(length(T2),1);
% for each T1 value, find the optimal base pays to maximize profit
for ii=1:length(T2)
    % find maximum profit by tuning tip rate, waiter base pay, and cook base pay
    [x,fval] = fminsearch(@(x)neg_profit(x,t0,tend,m1,m2,T1(ii),T2(ii),R,bw2,bc2,min_bw,min_bc,Nd,Nw,Nc,IC),initial_guess,optimset('TolFun',1e-4));
    % extract important parameters
    bw1(ii) = x(1); bc1(ii) = x(2); pmax(ii) = -fval;
end
end

function neg_P = neg_profit(x,t0,tend,m1,m2,T1,T2,R,bw2,bc2,min_bw,min_bc,Nd,Nw,Nc,IC)
% tune over waiter base pay x(1), and cook base pay x(2)
% impose huge penalty if parameters are out of bounds
if x(1)<min_bw || x(2)<min_bc
    neg_P = 10^10;
else
    % integrate system
    [~,Y] = integrate_tip_model(t0,tend,m1,m2,T1,T2,R,x(1),bw2,x(2),bc2,Nd,Nw,Nc,IC);

    % find equilibrium state
    stable_state = Y(end,:);

    % negative profit at equilibrium - ( + revenue - waiter pay - cook pay )
    neg_P = -stable_state(1)*m1 + stable_state(2)*x(1) + stable_state(3)*x(2);
end

end