% function returns best case equilibrium profit
function [pmax,T1,bw1,bc1] = max_profit3(t0,tend,m1,m2,T2,R,bw2,bc2,min_bw,min_bc,Nd,Nw,Nc,IC,initial_guess)
% 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,T2,R,bw2,bc2,min_bw,min_bc,Nd,Nw,Nc,IC),initial_guess,optimset('TolFun',1e-4));
% extract important parameters
T1 = x(1); bw1 = x(2); bc1 = x(3); pmax = -fval;

end

function neg_P = neg_profit(x,t0,tend,m1,m2,T2,R,bw2,bc2,min_bw,min_bc,Nd,Nw,Nc,IC)
% tune over tip rate x(1), waiter base pay x(2), and cook base pay x(3)
% impose huge penalty if parameters are out of bounds
if x(1)<0 || x(2)<min_bw || x(3)<min_bc
    neg_P = 10^10;
else

% integrate system
[~,Y] = integrate_tip_model(t0,tend,m1,m2,x(1),T2,R,x(2),bw2,x(3),bc2,Nd,Nw,Nc,IC);

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

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

end