function dydt=ODE_LHS(~,y,LHSmatrix,x,~) % PARAMETERS % Parameter_settings_LHS; R = LHSmatrix(x,1); Nc = LHSmatrix(x,2); Nd = LHSmatrix(x,3); m = LHSmatrix(x,4); T1 = LHSmatrix(x,5); T2 = LHSmatrix(x,6); bc1 = LHSmatrix(x,7); bc2 = LHSmatrix(x,8); bw1 = LHSmatrix(x,9); bw2 = LHSmatrix(x,10); C = Nc*bc1/(bc1+bc2); % assume C already at steady state for convenience Nw = 1; % y(1) is diners, y(2) is waiters, y(3) is cooks dydt = zeros(2,1); % value perceived by diners v1 = (y(2)+R*C)/(m*(1+T1)); v2 = (Nw-y(2)+R*(Nc-C))/(m*(1+T2)); % transition probabilities for diners Pd1 = v1/(v1+v2); Pd2 = v2/(v1+v2); % change in diners at our restaurant dydt(1) = (Nd-y(1))*Pd1 - y(1)*Pd2; % gratuity g1 = m*y(1)*T1/y(2); g2 = m*(Nd-y(1))*T2/(Nw-y(2)); % transition probabilities for waiters Pw1 = (bw1+g1)/(bw1+g1+bw2+g2); Pw2 = (bw2+g2)/(bw1+g1+bw2+g2); % change in waiters at our restaurant dydt(2) = (Nw-y(2))*Pw1 - y(2)*Pw2;