clc
clear

syms a b lambda theta_R Q c_S N n 

N = 5;
n = N - 1;
Q_R = sym ('q_R', [N 1]);
Q_R_non = Q_R;
Q_R_non(n) = 0;
W_R = sym ('w_R', [1 N]);
P_R = sym ('p_R', [N 1]);

% Price definition
for i = 1:N
    
    P_R(i) = a + Q_R(i) * (lambda - b) - lambda * Q;
    
end

p_M = a - b*Q + (b - lambda) * sum(Q_R);

p_bar_M = a - b*Q + Q_R(n)*(lambda - b) + sum(Q_R) * (b - lambda);

% Assumptions

assume( a > 0 );
assume( b > 0 );
assume( 0 <= lambda <= b);
assume( theta_R > 0 & theta_R < 1);
%theta_R=0;
assume( Q >= sum(Q_R) >= 0);
assume( c_S > 0);
assume( 0<= p_bar_M <= p_M );

for i=1:N
    
assume(Q_R(i) >= 0);
assume(0 < W_R(i) < P_R(i));

end

% Profit definition

PI_R = Q_R(n) * (P_R(n) - W_R(n));

PI_M = W_R*Q_R + (Q - sum(Q_R))* p_M - Q * c_S;

PI_bar_M = W_R * Q_R_non + (Q - sum(Q_R_non))* p_bar_M - Q * c_S;

% Nash Product

PHI = PI_R^(theta_R) * (PI_M - PI_bar_M)^(1 - theta_R);

PHI_diff = diff(PHI, Q_R(n));

x = solve(PHI_diff, Q_R(n))