%
% Note: This program does the calculation for both the case when b=0 (1D case,
%  the advertising costs are fixed)
%  and the case when the discount and advertising costs can vary.
%  set opt=1 to run the 1D case, opt=2 for the other
%
% Total profit from PC sales
%
 m = 700;          % manufacturing costs
 n = 10000;        % number sold w/o discount, extra ads
 p = 950;         % wholesale price
 r = .5;          % increase in sales per $100
 a = 50000;        %base advertising budget
 c = 200;          % increase in sales per $10,000 ads
% 1D - no increase in ads
 if opt == 1
 b = 0;
%
%First plot the function
%
 d1 = 0:.1:3;                     % discount range
 nn1= n.*(1+d1.*r) +c.*b;         % number sold due to discount
 q1 = (p-m)-d1.*100;              % profit/unit
 pp = nn1.*q1 - 50000 -b.*10000;  % total profit
 plot(pp,d1)
%
% symbolic differentiation, with maple
%
% symbolic variables 
 sr = sym('sr');                       % increase in sales per $100
 sd1 = sym('sd1');                     % discount (in $100) unit
 snn= n.*(1+sd1.*sr) +c.*b;            % number sold due to discount 
 sq1 = (p-m)-sd1.*100;                 % profit/unit
 spp = snn*sq1 - 50000 -b.*10000;      % total profit
%
% Differentiate: d(total profit)/d(Discount)
%
 dpd = diff(spp,sd1)                    
%
% Find optimal discount : set dpd=0 and solve for sd1
%
 [sd1] = solve('2500000*sr-2000000*sd1*sr-1000000','sd1')
%
% Now we have the optimal discount (sd1) as a function of the
% variable rate of increase of sales,  sr
%
% We now plot the discount for reasonable values of the rate r
% 
%
 rv = .41:.01:.6;                  % range of r
 dd1 = 1/4.*(5.*rv-2)./rv;         % optimal discount
 ddr = diff(sd1,sr) %(.5./sr.^2)   % d(optimal discount)/dr
 ddrv = subs(ddr,sr,rv)            % evaluate this derivative at rv
 sdr = ddrv.*rv./dd1;              % sensitivity function S(d,r)
 plot(rv,sdr)                      % 
 dpr = diff(spp,sr)                % d(total profit)/dr
 nnv= n.*(1+dd1.*rv) +c.*b;        % number sold, evaluated at rv 
 qv = (p-m)-dd1.*100;              % profit/unit, evaluated at rv
 ppv = nnv.*qv - 50000 -b.*10000;  % total profit, evaluated at rv
 dprv =-500000.*dd1.*(-5+2.*dd1)   % evaluate dp/dr at rv 
 spr = dprv.*rv./ppv;              % sensitivity function S(total profit, rv)
 plot(rv,spr)
%
%2D
%
else
clear pp;
%
% First plot the function
%
[d,b] = meshgrid(0:.025:3,0:.025:10);    % grid for discount and ad increase
 nn= n.*(1+d.*r) +c.*b;                  % number sold with discount /ads
 q= (p-m)-d.*100;                        % profit/unit
pp = nn.*q - 50000 -b.*10000;            % total profit
mesh(d,b,pp)                             % twod plot
%
% symbolic differentiation of the total profit
%
 sr = sym('sr');                        % rate of increase with discount 
 sd2 = sym('sd2');                      % discount
 sb = sym('sb');                        %  ad increase 
 snn= n.*(1+sd2.*sr) +c.*sb;            % number sold with discount/ads effect
 sq2 = (p-m)-sd2.*100;                  % number sold with discount
 spp = snn*sq2 - 50000 -sb.*10000;      % total profit
 dpd = diff(spp,sd2);                   % d(total profit)/d(discount)
 simplify(dpd)
 dpb = diff(spp,sb)                     % d(total profit)/d(b)
 simplify(dpb)
%
%  Solve system for dpd=0, dpb=0
%
 [sb,sd2]=solve('2500000*sr-2000000*sd2*sr-1000000-20000*sb=0','40000-20000*sd2 =0','sb,sd2')
end
%
% with lagrange multipliers
%
% First maximize on boundary b=5
%
 sd25 = sym('sd25');                    %  new discount variable
 lam5= sym('lam5');                     % new lambda
%
%  Solve system dpd=0 dpb = lam10 b=10 for discount, b, lambda
%
[b5,lam5,sd25] = solve('2500000*sr-2000000*sd25*sr-1000000-20000*b5=0', '40000-20000*sd25-lam5=0','b5=5','b5,lam5,sd25')
%
% Then  maximize on boudary b-0
% 
 sd20 = sym('sd20');
 lam0= sym('lam0');
%
%  Solve system dpd=0, dpb = lam10. b=0 for discount, b, lambda
%
[lam0,sd20] = solve('2500000*sr-2000000*sd20*sr-1000000-20000*0=0','40000-20000*sd20-lam0=0','lam0,sd20')