Stackerberg博弈：MATLAB代码实现与分析

%% 前言
%代码主题：stackerberg博弈，主要内容有：
%本代码基于win11系统下的MATLAB2016b运行，如代码运行失效，可能是软件版本不一致。
%请依次点击编辑器的【运行节】来运行代码。
%文中注释符：“%”、“%%”、“%{ %}”，您可根据需要自行删除注释。
%原论文中的希腊字母均用其发音表示，具体可查看百度，其他字母直接观察即可得知其原文中的对应变量名。
%论文作者与标题：李重莲等，《双向公平关切下双渠道供应链的线上线下融合契约设计》。
%代码作者：b站UP主【谭小飞導仕】（原名：谭小飞同学），如有疑问或是其他问题，可在b站联系。
%代码完成时间：2023年3月5日。
%% 3.1 双渠道供应链的集中决策模型
%{
思路：
1将制造商、线上渠道商和线下渠道商视为同一企业的不同部门，先求出该企业的总期望净收益PIc，
并分别列出线上渠道商的需求函数do与线下渠道商的需求函数dt；
2其次，证明收益函数PIc分别对线上渠道商零售价格po和线下渠道商零售价格pt的海塞矩阵均为负，
于是，存在唯一的最优解。那么，再以收益函数PIc分别对po、pt求一阶导并为零，得出最优解po*、pt*；
3再次，将最优零售价po*、pt*分别代入需求函数do、dt中即可得知最优产品需求量do*、dt*；
4最后，将po*、pt*、do*、dt*全部代入PIc中可得最大化的期望收益PIc*。
注意：3.1节无需求出w，因为将上述直接带入PIc后，w的乘数项为零，因此无需求制造商批发价格w。
%}

clc,clear;%清空命令行窗口，清空工作区内所有变量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%第一步：列出系统总收益函数以及渠道商的需求函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms PIc w c Q beta pt po co k ct dt do; % 定义字符型变量
PIc=(w-c)*(Q-beta*(pt+po))+(po-w-co)*((1-k)*Q-beta*po)+(pt-w-ct)*(k*Q-beta*pt);% 系统的总的净收益函数
dt=k*Q-beta*pt;%线下渠道商实际需求函数
do=(1-k)*Q-beta*po;%线上渠道商实际需函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%第二步：求出海塞矩阵，并求出一阶导对应的最优解%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%先求出海塞矩阵是否为负定。
H=hessian(PIc,[po;pt]);%求PIc对变量po、pt的二阶偏导，即海塞矩阵，变量间以分号隔开，矩阵记为H。
fprintf('海塞矩阵H=\n');%在命令行窗口显示海塞矩阵H的结果。
disp(H);%在命令行窗口显示矩阵H。
fprintf('海塞矩阵H的latex代码为：\n');%在命令行窗口显示矩阵H对应的latex代码。
disp(latex(H)); %在命令行窗口显示矩阵H的latex代码，可直接粘贴至latex或mathtype中。
H = subs(H, beta,1); %MATLAB只能比较数值型而不能比较字符型变量的大小，因论文设定beta>0，故这里随意地赋值为1，以方便后面矩阵判定。
d = eig(H); % 计算矩阵特征值
if all(d > 0)
    fprintf('海塞矩阵H正定.\n');%如果海塞矩阵特征值均为正，则为正定矩阵。
elseif all(d >= 0)
    fprintf('海塞矩阵H半正定.\n');%如果海塞矩阵特征值均为非负，则为半正定矩阵。
else
    fprintf('海塞矩阵H负定。\n');%如果海塞矩阵特征值均为负，则为负定矩阵，意味着函数PIc关于变量po、pt存在唯一一组最大值解。
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%在海塞矩阵判定为负的基础上，再求出最优解po*、pt*。
%%%%%%%%%%%%%% 1 求出po*、pt*；
equ1 = diff(PIc, po);%PIc对po求一阶偏导，将求导结果存储在暂元equ1中。
poc = solve(equ1, po);%上式equ1为零时po对应的解，记为最优解poc
equ2 = diff(PIc, pt);%PIc对pt求一阶偏导，将求导结果存储在暂元equ2中。
ptc = solve(equ2, pt);%上式equ2为零时pt对应的解，记为最优解ptc
%%%%%%%%%%%%%% 2 显示po*；
fprintf(2,'*1最优解poc为*：');%在命令行窗口以红色字体突出显示。
fprintf('poc=');
disp(simplify(poc));
fprintf('poc的数学表达式为：\n');
pretty(simplify(poc));%将最优解poc化简，并以手写形式显示
fprintf('poc的latex代码为：\n');
disp(latex(simplify(poc)));
%%%%%%%%%%%%%% 3 显示pt*；
fprintf(2,'*2最优解ptc为*：');
fprintf('ptc=');
disp(simplify(ptc));
fprintf('ptc的数学表达式为：\n');
pretty(simplify(ptc));
fprintf('ptc的latex代码为：\n');
disp(latex(simplify(ptc)));
%%%%%%%%%%%%%% 4 将po*、pt*分别代入需求函数do、dt中；
dtc = subs(dt, pt, ptc);
doc = subs(do, po, poc);
%%%%%%%%%%%%%% 5 显示dt*；
fprintf(2,'*3最优解dtc为*：');
fprintf('dtc=');
disp(simplify(dtc));
fprintf('dtc化简后的数学表达式为：\n');
pretty(simplify(dtc));
fprintf('dtc的latex代码为：\n');
disp(latex(simplify(dtc)));
%%%%%%%%%%%%%% 6 显示do*；
fprintf(2,'*4最优解doc为*：');
fprintf('doc=');
disp(simplify(doc));
fprintf('doc化简后的数学表达式为：\n');
pretty(simplify(doc));
fprintf('doc的latex代码为：\n');
disp(latex(simplify(doc)));
%%%%%%%%%%%%%% 7 将po*、pt*、dt*、do*代入收益函数PI中；
PIc = subs(PIc, {pt,po,do,dt},{ptc,poc,doc,dtc});
%%%%%%%%%%%%%% 8 显示PI*；
fprintf(2,'*5最优解PIc为*：');
fprintf('PIc=');
disp(simplify(PIc));
fprintf('供应链的总利润PIc为：\n');
pretty(simplify(PIc));
fprintf('PIc的latex代码为：\n');
disp(latex(simplify(PIc)));
%% 3.2 双渠道供应链的分散决策模型
%{
思路：
1.在分散决策模型中，制造商、线下渠道商和线上渠道商是不同的企业，以制造商为主、以渠道商为从，
先列出三者的期望收益函数。
2.分别求出线上、线下制造商收益函数对各自零售价格的海塞矩阵，然后再求出一阶导的最优价格。
其实，在该模型中，并不需要求出海塞矩阵，因为海塞矩阵判定是用于多元函数如f(x,y)分别对x、y求二阶偏导，
本模型中的制造商收益函数是一元二次函数，并不需要用到海塞矩阵，但为了与原文保持一致，本人仍提供相应代码。
3.将求出的渠道商的最优零售价格代入制造商的收益函数中，并求制造商收益函数对其批发价格w的一阶导，得到最优批发价。
4.用字母G替换字母w，虽然本人觉得多此一举，但仍照做。
5.将最优零售价分别代入需求函数中即可得知最优产品需求量；
6.最后，将以上所有最优零售价、最优批发价以及最优需求量全部代入相应的收益函数中。
%}
clc,clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1 列出所有收益函数与渠道商需求函数。
syms w c Q beta pt po co k ct pt; %换了模型后，本人清除了前面所有结果，这里重新定义。
PImw=(w-c)*(Q-beta*(pt+po));
PIow=(po-w-co)*((1-k)*Q-beta*po);
PItw=(pt-w-ct)*(k*Q-beta*pt);
dt=k*Q-beta*pt;
do=(1-k)*Q-beta*po;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2 判定两个海塞矩阵是否为负定。
H_PIow=hessian(PIow,po);%求PIow的二阶海塞矩阵
fprintf('海塞矩阵H_PIow为=\n');
disp(H_PIow);
H_PItw=hessian(PItw,pt);%求PImw的二阶海塞矩阵
fprintf('海塞矩阵H_PItw为=\n');
disp(H_PItw);
H_PIow = subs(H_PIow, beta,1);
d = eig(H_PIow); % 计算矩阵特征值
if all(d > 0)
    fprintf('海塞矩阵H_PIow正定.\n');
elseif all(d >= 0)
    fprintf('海塞矩阵H_PIow半正定.\n');
else
    fprintf('海塞矩阵H_PIow负定。\n');
end
H_PItw = subs(H_PItw, beta,1);
d = eig(H_PItw); % 计算矩阵特征值
if all(d > 0)
    fprintf('海塞矩阵H_PItw正定.\n');
elseif all(d >= 0)
    fprintf('海塞矩阵H_PItw半正定.\n');
else
    fprintf('海塞矩阵H_PItw负定。\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 求出最优解po*、pt*。
equ1 = diff(PIow, po);
p_o = solve(equ1, po);
equ2 = diff(PItw, pt);
p_t = solve(equ2, pt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4 分别列出最优解po*、pt*。
fprintf(2,'*1最优解p_o为*：');
fprintf('p_o=');
disp(simplify(p_o));
fprintf(2,'*2最优解p_t为*：');
fprintf('p_t=');
disp(simplify(p_t));
fprintf('最优解p_o、p_t分别为：\n');
pretty(simplify(p_o)),pretty(simplify(p_t));
fprintf('p_o的latex代码为：\n');
disp(latex(simplify(p_o)));
fprintf('p_t的latex代码为：\n');
disp(latex(simplify(p_t)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5 将最优解po*、pt*代入制造商收益函数PImw中，并求PImw对变量w的一阶导的最优解w*。
PImw = subs(PImw, {po,pt},{p_o,p_t});
equ = diff(PImw, w);
ww = solve(equ, w);
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 6 列出最优解w*。
fprintf(2,'*3最优解ww为*：');
fprintf('ww=');
disp(simplify(ww));
fprintf('最优解ww的解为：\n');
pretty(simplify(ww));
fprintf('ww的latex代码为：\n');
disp(latex(simplify(ww)));
%%%%%%%%%%%%%% 7 将po*、pt*、w*分别代入收益函数PI、需求函数d中；
PImw = subs(PImw, {pt,po,w},{p_t,p_o,ww});
PIow = subs(PIow, {po,w},{p_o,ww});
PItw = subs(PItw, {pt,w},{p_t,ww});
dtw = subs(dt, pt,p_t);
dow = subs(do, po,p_o);
%%%%%%%%%%%%%% 8 显示PI*、d*；
fprintf(2,'*4最优解PImw为*：');
fprintf('PImw=');
disp(simplify(PImw));
fprintf('制造商收益PImw为：\n');
pretty(simplify(PImw));
fprintf('PImw的latex代码为：\n');
disp(latex(simplify(PImw)));

fprintf(2,'*5最优解PIow为*：');
fprintf('PIow=');
disp(simplify(PIow));
fprintf('制造商收益PIow为：\n');
pretty(simplify(PIow));
fprintf('PIow的latex代码为：\n');
disp(latex(simplify(PIow)));

fprintf(2,'*6最优解PItw为*：');
fprintf('PItw=');
disp(simplify(PItw));
fprintf('制造商收益PItw为：\n');
pretty(simplify(PItw));
fprintf('PItw的latex代码为：\n');
disp(latex(simplify(PItw)));

fprintf(2,'*7最优解dtw为*：');
fprintf('dtw=');
disp(simplify(dtw));
fprintf('线下渠道商销售量dtw为：\n');
pretty(simplify(dtw));
fprintf('dtw的latex代码为：\n');
disp(latex(simplify(dtw)));

fprintf(2,'*8最优解dow为*：');
fprintf('dow=');
disp(simplify(dow));
fprintf('线上渠道商销售量dow为：\n');
pretty(simplify(dow));
fprintf('dow的latex代码为：\n');
disp(latex(simplify(dow)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 9 求出供应链的总收益。
PIw=PImw+PIow+PItw;
fprintf(2,'*9最优解PIw为*：');
fprintf('PIw=');
disp(simplify(PIw));
fprintf('供应链总收益PIw为：\n');
pretty(simplify(PIw));
fprintf('PIw的latex代码为：\n');
disp(latex(simplify(PIw)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 10 用字母G替换字母w，这一步可以省略，因为原文后面赋值时没有考虑G。
fprintf(2,'-----------------------------------------------------我是分割线-----------------------------------------------------\n');
fprintf(2,'-----------------------------------------------------我是分割线-----------------------------------------------------\n');
fprintf(2,'-----------------------------------------------------我是分割线-----------------------------------------------------\n');
syms Gw;
pow=subs(p_o,w,Gw);
ptw=subs(p_t,w,Gw);
dow=subs(do,po,pow);
dtw=subs(dt,pt,ptw);
fprintf('最优零售价pow、ptw分别为：\n');
pretty(simplify(pow)),pretty(simplify(ptw));
fprintf('pow、ptw的latex代码为：\n');
disp(latex(simplify(pow))),disp(latex(simplify(ptw)));
fprintf('dow、dtw分别为：\n');
pretty(simplify(dow)),pretty(simplify(dtw));
fprintf('dow、dtw的latex代码为：\n');
disp(latex(simplify(dow))),disp(latex(simplify(dtw)));
%% 3.2 双渠道供应链的公平关切下的分散决策模型
%{
思路：
本小节与上一小节基本类似，只是多引入了两个公平关切系数gamma和delta，收益函数不再用PI而是U表示。
%}
clc,clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1 列出期望收益和期望效用函数、需求函数。
syms w c Q beta pt po co k ct pt delta gamma;
PIm=(w-c)*(Q-beta*(pt+po));
PIo=(po-w-co)*((1-k)*Q-beta*po);
PIt=(pt-w-ct)*(k*Q-beta*pt);
Uon=PIo-gamma*(PIt-PIo)-delta*(PIm-PIo);
Utn=PIt-gamma*(PIo-PIt)-delta*(PIm-PIt);
Umn=PIm-delta*(PIo-PIm)-delta*(PIt-PIm);
dt=k*Q-beta*pt;
do=(1-k)*Q-beta*po;
pretty(simplify(Umn));
pretty(simplify(Utn));
pretty(simplify(Uon));
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2 求出并显示最优解po*、pt*。
equ1 = diff(Uon, po);
pon = solve(equ1, po);
equ2 = diff(Utn, pt);
ptn = solve(equ2, pt);
fprintf(2,'*1最优解pon为*：');
fprintf('pon=');
disp(simplify(pon));
fprintf(2,'*2最优解ptn为*：');
fprintf('ptn=');
disp(simplify(ptn));
fprintf('最优解pon、ptn分别为：\n');
pretty(simplify(pon)),pretty(simplify(ptn));
fprintf('最优解pon、ptn的latex代码分别为：\n');
disp(latex(simplify(pon))),disp(latex(simplify(ptn)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 求出并显示最优解w*。
Umn = subs(Umn, {po,pt},{pon,ptn});
equ = diff(Umn, w);
wn = solve(equ, w);
fprintf(2,'*3最优解wn为*：');
fprintf('wn=');
disp(simplify(wn));
fprintf('最优解wn的解为：\n');
pretty(simplify(wn));
fprintf('wn的latex代码为：\n');
disp(latex(simplify(wn)));
%%%%%%%%%%%%%% 4 将po*、pt*、w*分别代入收益函数U、需求函数d中；
Umn = subs(Umn, {pt,po,w},{ptn,pon,wn});
Uon = subs(Uon, {pt,po,w},{ptn,pon,wn});
Utn = subs(Utn, {pt,po,w},{ptn,pon,wn});
dtn = subs(dt, pt,ptn);
don = subs(do, po,pon);
%%%%%%%%%%%%%% 5 显示U*、d*；
fprintf(2,'*4最优解Umn为*：');
fprintf('Umn=');
disp(simplify(Umn));
fprintf('制造商收益Umn为：\n');
pretty(simplify(Umn));
fprintf('Umn的latex代码为：\n');
disp(latex(simplify(Umn)));
fprintf(2,'*5最优解Uon为*：');
fprintf('Uon=');
disp(simplify(Uon));
fprintf('制造商收益Uon为：\n');
pretty(simplify(Uon));
fprintf('Uon的latex代码为：\n');
disp(latex(simplify(Uon)));
fprintf(2,'*6最优解Utn为*：');
fprintf('Utn=');
disp(simplify(Utn));
fprintf('制造商收益Utn为：\n');
pretty(simplify(Utn));
fprintf('Utn的latex代码为：\n');
disp(latex(simplify(Utn)));
fprintf(2,'*7最优解dtn为*：');
fprintf('dtn=');
disp(simplify(dtn));
fprintf('线下渠道商销售量dtn为：\n');
pretty(simplify(dtn));
fprintf('dtn的latex代码为：\n');
disp(latex(simplify(dtn)));
fprintf(2,'*8最优解don为*：');
fprintf('don=');
disp(simplify(don));
fprintf('线上渠道商销售量don为：\n');
pretty(simplify(don));
fprintf('don的latex代码为：\n');
disp(latex(simplify(don)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 6 求出供应链的总收益。
Un=Umn+Uon+Utn;
fprintf(2,'*9最优解Un为*：');
fprintf('Un=');
disp(simplify(Un));
fprintf('供应链总收益Un为：\n');
pretty(simplify(Un));
fprintf('Un的latex代码为：\n');
disp(latex(simplify(Un)));
%% 3.3 双向公平关切行为对双渠道供应链决策的影响分析
%{
由于论文中的本小节部分是比较字符型变量之间的大小，且涉及数学推导过程，因此冀希望于MATLAB是不现实的，故本节无代码，特此说明。
%}
%% 4.1 契约模型的建立
%{
思路：
在公平关切下的分散决策模型基础上引入参数theta、phi。
不同于3.2节可以通过函数代入来减少输入时间，但本小节的收益函数需要逐个输入。
%}
clc,clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1 列出所有收益函数与渠道商需求函数。
syms theta phi delta w c Q beta pt po ct co k gamma;
PImx=(1+2*delta)*(theta*w-c)*(Q-beta*(pt+po))-delta*(pt-theta*w-(1-phi)*ct)*(k*Q-beta*pt)-delta*(po-theta*w-co-phi*ct)*((1-k)*Q-beta*po);
PIox=(1+gamma+delta)*(po-theta*w-co-phi*ct)*((1-k)*Q-beta*po)-gamma*(pt-theta*w-(1-phi)*ct)*(k*Q-beta*pt)-delta*(theta*w-c)*(Q-beta*(pt+po));
PItx=(1+gamma+delta)*(pt-theta*w-(1-phi)*ct)*(k*Q-beta*pt)-gamma*(po-theta*w-co-phi*ct)*((1-k)*Q-beta*po)-delta*(theta*w-c)*(Q-beta*(pt+po));
PIx=PImx+PIox+PItx;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2 求出并显示最优解po*、pt*。
equ1 = diff(PIox, po);
pox = solve(equ1, po);
equ2 = diff(PItx, pt);
ptx = solve(equ2, pt);
fprintf(2,'*1最优解pox为*：');
fprintf('pox=');
disp(simplify(pox));
fprintf(2,'*2最优解ptx为*：');
fprintf('ptx=');
disp(simplify(ptx));
fprintf('最优解pox、ptx分别为：\n');
pretty(simplify(pox)),pretty(simplify(ptx));
fprintf('最优解pox、ptx的latex代码分别为：\n');
disp(latex(simplify(pox))),disp(latex(simplify(ptx)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 求出并显示最优解w*。
%根据论文解释，w*的求解步骤是：先求出满足等式poc=pox、ptc=ptx的w1、w2，再令w*=max{w1,w2}。
poc=(Q - Q*k + beta*c + beta*co)/(2*beta);
ptc=(Q*k + beta*c + beta*ct)/(2*beta);
pox=-(beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1));
ptx=(beta*(theta*w - ct*(phi - 1))*(delta + gamma + 1) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1));
w1 = solve(poc-pox, w);
w2 = solve(ptc-ptx, w);
fprintf(2,'w1=\n');pretty(simplify(w1));
fprintf(2,'w2=\n');pretty(simplify(w2));
%由于MATLAB不能比较字符型变量的大小，而赋值比较又较为麻烦，因此在命令行窗口直接观察即可得知w2>w1，则wx=w2.
wx=w2;
fprintf(2,'*3最优解wx为*：');
fprintf('wx=');
disp(simplify(wx));
%%%%%%%%%%%%%% 4 将po*、pt*、w*分别代入收益函数PI中；
PImx = subs(PImx, {pt,po,w},{ptx,pox,wx});
PIox = subs(PIox, {pt,po,w},{ptx,pox,wx});
PItx = subs(PItx, {pt,po,w},{ptx,pox,wx});
%%%%%%%%%%%%%% 5 显示PI*；
fprintf(2,'*4最优解PImx为*：');
fprintf('PImx=');
disp(simplify(PImx));
fprintf('制造商收益PImx为：\n');
pretty(simplify(PImx));
fprintf('PImx的latex代码为：\n');
disp(latex(simplify(PImx)));
fprintf(2,'*5最优解PIox为*：');
fprintf('PIox=');
disp(simplify(PIox));
fprintf('线上制造商收益PIox为：\n');
pretty(simplify(PIox));
fprintf('PIox的latex代码为：\n');
disp(latex(simplify(PIox)));
fprintf(2,'*6最优解PItx为*：');
fprintf('PItx=');
disp(simplify(PItx));
fprintf('线下制造商收益PItx为：\n');
pretty(simplify(PItx));
fprintf('PItx的latex代码为：\n');
disp(latex(simplify(PItx)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 6 求出供应链PIx的总收益。
PIx=PImx+PIox+PItx;
fprintf(2,'*7最优解PIx为*：');
fprintf('PIx=');
disp(simplify(PIx));
fprintf('供应链总收益PIx为：\n');
pretty(simplify(PIx));
fprintf('PIx的latex代码为：\n');
disp(latex(simplify(PIx)));
%% 4.2 双渠道供应链契约协调的效果分析
%{
本小节属于数学推导，同样不提供MATLAB代码。原因同上。
%}
%% 5.1 三种双渠道供应链决策模型的对比分析
clc,clear;
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8;%参数赋值
%%%%%%%%%%%%%%%%%%%%%%%%1集中决策的各项参数值%%%%%%%%%%%%%%%%%%%%%%%%
poc=(Q - Q*k + beta*c + beta*co)/(2*beta);
ptc=(Q*k + beta*c + beta*ct)/(2*beta);
dtc=(Q*k)/2 - (beta*c)/2 - (beta*ct)/2;
doc=Q/2 - (Q*k)/2 - (beta*c)/2 - (beta*co)/2;
PIc=(2*Q^2*k^2 - 2*Q^2*k + Q^2 - 2*Q*beta*c + 2*Q*beta*co*k - 2*Q*beta*co - 2*Q*beta*ct*k + 2*beta^2*c^2 + 2*beta^2*c*co + 2*beta^2*c*ct + beta^2*co^2 + beta^2*ct^2)/(4*beta);
fprintf('-----------------------------------------------------我是一号分割线-----------------------------------------------------\n');
fprintf(2,'wc=-');
poc=roundn(poc,-2);fprintf(2,['\npoc=',num2str(poc)]);
ptc=roundn(ptc,-2);fprintf(2,['\nptc=',num2str(ptc)]);
doc=roundn(doc,-2);fprintf(2,['\ndoc=',num2str(doc)]);
dtc=roundn(dtc,-2);fprintf(2,['\ndtc=',num2str(dtc)]);%经过笔者手算，论文中是错误的结果。
fprintf(2,'\nPImc=-');
fprintf(2,'\nPIoc=-');
fprintf(2,'\nPItc=-');
PIc=roundn(PIc,-2);fprintf(2,['\nPIc=',num2str(PIc)]);%虽然数字与论文接近，且笔者未经过手算，但笔者仍认为是论文计算错误。
disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%2分散决策的各项参数值%%%%%%%%%%%%%%%%%%%%%%%%
ww=(Q + 2*beta*c - beta*co - beta*ct)/(4*beta);
w=ww;
p_o=(Q - Q*k + beta*co + beta*w)/(2*beta);
p_t=(Q*k + beta*(ct + w))/(2*beta);

PImw=(2*beta*c - Q + beta*co + beta*ct)^2/(16*beta);
PIow=-(((Q*k)/2 - Q/2 + (beta*co)/2 + (beta*w)/2)*(Q - 2*Q*k - 2*beta*c - beta*co + beta*ct + 2*beta*w))/(4*beta);
PItw=(((beta*ct)/2 - (Q*k)/2 + (beta*w)/2)*(Q - 2*Q*k + 2*beta*c - beta*co + beta*ct - 2*beta*w))/(4*beta);
dtw=(Q*k)/2 - (beta*(ct + w))/2;
dow=Q/2 - (Q*k)/2 - (beta*co)/2 - (beta*w)/2;
PIw=(8*Q^2*k^2 - 8*Q^2*k + 3*Q^2 - 8*Q*beta*c + 8*Q*beta*co*k - 6*Q*beta*co - 8*Q*beta*ct*k + 2*Q*beta*ct + 4*Q*beta*w + 4*beta^2*c^2 + 8*beta^2*c*co + 8*beta^2*c*ct + 8*beta^2*c*w + 3*beta^2*co^2 - 2*beta^2*co*ct - 4*beta^2*co*w + 3*beta^2*ct^2 - 4*beta^2*ct*w - 8*beta^2*w^2)/(16*beta);
fprintf('-----------------------------------------------------我是二号分割线-----------------------------------------------------\n');
ww=roundn(ww,-2);fprintf(2,['ww=',num2str(ww)]);
p_o=roundn(p_o,-2);fprintf(2,['\np_o=',num2str(p_o)]);
p_t=roundn(p_t,-2);fprintf(2,['\np_t=',num2str(p_t)]);
dow=roundn(dow,-2);fprintf(2,['\ndow=',num2str(dow)]);
dtw=roundn(dtw,-2);fprintf(2,['\ndtw=',num2str(dtw)]);
PImw=roundn(PImw,-2);fprintf(2,['\nPImw=',num2str(PImw)]);
PIow=roundn(PIow,-2);fprintf(2,['\nPIow=',num2str(PIow)]);
PItw=roundn(PItw,-2);fprintf(2,['\nPItw=',num2str(PItw)]);
PIw=roundn(PIw,-2);fprintf(2,['\nPIw=',num2str(PIw)]);
disp(' ');
%%%%%%%%%%%%%%%%%%%%%%%%3关切下分散决策的各项参数值%%%%%%%%%%%%%%%%%%%%%%%%
wn=(Q + 5*Q*delta + 2*Q*gamma + 2*beta*c - beta*co - beta*ct + 7*Q*delta^2 + 3*Q*delta^3 + Q*gamma^2 + 8*Q*delta*gamma + 12*beta*c*delta - 5*beta*co*delta - 5*beta*ct*delta + 4*beta*c*gamma - 2*beta*co*gamma - 2*beta*ct*gamma + 3*Q*delta*gamma^2 + 6*Q*delta^2*gamma + 22*beta*c*delta^2 + 10*beta*c*delta^3 - 7*beta*co*delta^2 - 3*beta*co*delta^3 - 7*beta*ct*delta^2 - 3*beta*ct*delta^3 + 2*beta*c*gamma^2 - beta*co*gamma^2 - beta*ct*gamma^2 + 4*beta*c*delta*gamma^2 + 16*beta*c*delta^2*gamma - 3*beta*co*delta*gamma^2 - 6*beta*co*delta^2*gamma - 3*beta*ct*delta*gamma^2 - 6*beta*ct*delta^2*gamma + 16*beta*c*delta*gamma - 8*beta*co*delta*gamma - 8*beta*ct*delta*gamma)/(2*beta*(8*delta^3 + 14*delta^2*gamma + 18*delta^2 + 5*delta*gamma^2 + 16*delta*gamma + 11*delta + 2*gamma^2 + 4*gamma + 2));
w=wn;
pon=-(Q*(k - 1) + gamma*(Q*(k - 1) - beta*(co + w)) + delta*(beta*(c - w) + Q*(k - 1) - beta*(co + w)) - beta*(co + w))/(2*beta + 2*beta*delta + 2*beta*gamma);
ptn=(Q*k + gamma*(Q*k + beta*(ct + w)) + delta*(Q*k - beta*(c - w) + beta*(ct + w)) + beta*(ct + w))/(2*beta + 2*beta*delta + 2*beta*gamma);
don=(beta*(Q*(k - 1) + gamma*(Q*(k - 1) - beta*(co + w)) + delta*(beta*(c - w) + Q*(k - 1) - beta*(co + w)) - beta*(co + w)))/(2*beta + 2*beta*delta + 2*beta*gamma) - Q*(k - 1);
dtn=Q*k - (beta*(Q*k + gamma*(Q*k + beta*(ct + w)) + delta*(Q*k - beta*(c - w) + beta*(ct + w)) + beta*(ct + w)))/(2*beta + 2*beta*delta + 2*beta*gamma);
Umn=(32*Q^2*delta^4*k - 32*Q^2*delta^4*k^2 - 7*Q^2*delta^4 - 56*Q^2*delta^3*gamma*k^2 + 56*Q^2*delta^3*gamma*k - 10*Q^2*delta^3*gamma - 72*Q^2*delta^3*k^2 + 72*Q^2*delta^3*k - 12*Q^2*delta^3 - 20*Q^2*delta^2*gamma^2*k^2 + 20*Q^2*delta^2*gamma^2*k - Q^2*delta^2*gamma^2 - 64*Q^2*delta^2*gamma*k^2 + 64*Q^2*delta^2*gamma*k - 2*Q^2*delta^2*gamma - 44*Q^2*delta^2*k^2 + 44*Q^2*delta^2*k - 8*Q^2*delta*gamma^2*k^2 + 8*Q^2*delta*gamma^2*k + 2*Q^2*delta*gamma^2 - 16*Q^2*delta*gamma*k^2 + 16*Q^2*delta*gamma*k + 6*Q^2*delta*gamma - 8*Q^2*delta*k^2 + 8*Q^2*delta*k + 4*Q^2*delta + Q^2*gamma^2 + 2*Q^2*gamma + Q^2 - 4*Q*beta*c*delta^4 - 16*Q*beta*c*delta^3*gamma - 24*Q*beta*c*delta^3 - 16*Q*beta*c*delta^2*gamma^2 - 56*Q*beta*c*delta^2*gamma - 44*Q*beta*c*delta^2 - 16*Q*beta*c*delta*gamma^2 - 40*Q*beta*c*delta*gamma - 24*Q*beta*c*delta - 4*Q*beta*c*gamma^2 - 8*Q*beta*c*gamma - 4*Q*beta*c - 32*Q*beta*co*delta^4*k + 14*Q*beta*co*delta^4 - 56*Q*beta*co*delta^3*gamma*k + 20*Q*beta*co*delta^3*gamma - 72*Q*beta*co*delta^3*k + 24*Q*beta*co*delta^3 - 20*Q*beta*co*delta^2*gamma^2*k + 2*Q*beta*co*delta^2*gamma^2 - 64*Q*beta*co*delta^2*gamma*k + 4*Q*beta*co*delta^2*gamma - 44*Q*beta*co*delta^2*k - 8*Q*beta*co*delta*gamma^2*k - 4*Q*beta*co*delta*gamma^2 - 16*Q*beta*co*delta*gamma*k - 12*Q*beta*co*delta*gamma - 8*Q*beta*co*delta*k - 8*Q*beta*co*delta - 2*Q*beta*co*gamma^2 - 4*Q*beta*co*gamma - 2*Q*beta*co + 32*Q*beta*ct*delta^4*k - 18*Q*beta*ct*delta^4 + 56*Q*beta*ct*delta^3*gamma*k - 36*Q*beta*ct*delta^3*gamma + 72*Q*beta*ct*delta^3*k - 48*Q*beta*ct*delta^3 + 20*Q*beta*ct*delta^2*gamma^2*k - 18*Q*beta*ct*delta^2*gamma^2 + 64*Q*beta*ct*delta^2*gamma*k - 60*Q*beta*ct*delta^2*gamma + 44*Q*beta*ct*delta^2*k - 44*Q*beta*ct*delta^2 + 8*Q*beta*ct*delta*gamma^2*k - 12*Q*beta*ct*delta*gamma^2 + 16*Q*beta*ct*delta*gamma*k - 28*Q*beta*ct*delta*gamma + 8*Q*beta*ct*delta*k - 16*Q*beta*ct*delta - 2*Q*beta*ct*gamma^2 - 4*Q*beta*ct*gamma - 2*Q*beta*ct + 4*beta^2*c^2*delta^4 + 16*beta^2*c^2*delta^3*gamma + 24*beta^2*c^2*delta^3 + 16*beta^2*c^2*delta^2*gamma^2 + 56*beta^2*c^2*delta^2*gamma + 44*beta^2*c^2*delta^2 + 16*beta^2*c^2*delta*gamma^2 + 40*beta^2*c^2*delta*gamma + 24*beta^2*c^2*delta + 4*beta^2*c^2*gamma^2 + 8*beta^2*c^2*gamma + 4*beta^2*c^2 + 4*beta^2*c*co*delta^4 + 16*beta^2*c*co*delta^3*gamma + 24*beta^2*c*co*delta^3 + 16*beta^2*c*co*delta^2*gamma^2 + 56*beta^2*c*co*delta^2*gamma + 44*beta^2*c*co*delta^2 + 16*beta^2*c*co*delta*gamma^2 + 40*beta^2*c*co*delta*gamma + 24*beta^2*c*co*delta + 4*beta^2*c*co*gamma^2 + 8*beta^2*c*co*gamma + 4*beta^2*c*co + 4*beta^2*c*ct*delta^4 + 16*beta^2*c*ct*delta^3*gamma + 24*beta^2*c*ct*delta^3 + 16*beta^2*c*ct*delta^2*gamma^2 + 56*beta^2*c*ct*delta^2*gamma + 44*beta^2*c*ct*delta^2 + 16*beta^2*c*ct*delta*gamma^2 + 40*beta^2*c*ct*delta*gamma + 24*beta^2*c*ct*delta + 4*beta^2*c*ct*gamma^2 + 8*beta^2*c*ct*gamma + 4*beta^2*c*ct - 7*beta^2*co^2*delta^4 - 10*beta^2*co^2*delta^3*gamma - 12*beta^2*co^2*delta^3 - beta^2*co^2*delta^2*gamma^2 - 2*beta^2*co^2*delta^2*gamma + 2*beta^2*co^2*delta*gamma^2 + 6*beta^2*co^2*delta*gamma + 4*beta^2*co^2*delta + beta^2*co^2*gamma^2 + 2*beta^2*co^2*gamma + beta^2*co^2 + 18*beta^2*co*ct*delta^4 + 36*beta^2*co*ct*delta^3*gamma + 48*beta^2*co*ct*delta^3 + 18*beta^2*co*ct*delta^2*gamma^2 + 60*beta^2*co*ct*delta^2*gamma + 44*beta^2*co*ct*delta^2 + 12*beta^2*co*ct*delta*gamma^2 + 28*beta^2*co*ct*delta*gamma + 16*beta^2*co*ct*delta + 2*beta^2*co*ct*gamma^2 + 4*beta^2*co*ct*gamma + 2*beta^2*co*ct - 7*beta^2*ct^2*delta^4 - 10*beta^2*ct^2*delta^3*gamma - 12*beta^2*ct^2*delta^3 - beta^2*ct^2*delta^2*gamma^2 - 2*beta^2*ct^2*delta^2*gamma + 2*beta^2*ct^2*delta*gamma^2 + 6*beta^2*ct^2*delta*gamma + 4*beta^2*ct^2*delta + beta^2*ct^2*gamma^2 + 2*beta^2*ct^2*gamma + beta^2*ct^2)/(8*beta*(8*delta^3 + 14*delta^2*gamma + 18*delta^2 + 5*delta*gamma^2 + 16*delta*gamma + 11*delta + 2*gamma^2 + 4*gamma + 2));
Uon=-(3*Q^2*k - 5*Q^2*gamma - 8*Q^2*delta - Q^2 - 25*Q^2*delta^2 - 39*Q^2*delta^3 - 32*Q^2*delta^4 - 13*Q^2*delta^5 - 2*Q^2*delta^6 - 10*Q^2*gamma^2 - 10*Q^2*gamma^3 - 5*Q^2*gamma^4 - Q^2*gamma^5 - 2*Q^2*k^2 - beta^2*co^2 + 2*beta^2*w^2 - 56*Q^2*delta*gamma^2 - 86*Q^2*delta^2*gamma - 44*Q^2*delta*gamma^3 - 101*Q^2*delta^3*gamma - 16*Q^2*delta*gamma^4 - 55*Q^2*delta^4*gamma - 2*Q^2*delta*gamma^5 - 11*Q^2*delta^5*gamma - 17*Q^2*delta*k^2 + 89*Q^2*delta^2*k + 154*Q^2*delta^3*k + 143*Q^2*delta^4*k + 68*Q^2*delta^5*k + 13*Q^2*delta^6*k - 8*Q^2*gamma*k^2 + 26*Q^2*gamma^2*k + 24*Q^2*gamma^3*k + 11*Q^2*gamma^4*k + 2*Q^2*gamma^5*k + 2*beta^2*c^2*delta - 8*beta^2*co^2*delta + beta^2*ct^2*delta - 5*beta^2*co^2*gamma + beta^2*ct^2*gamma + 21*beta^2*delta*w^2 + 8*beta^2*gamma*w^2 + 2*Q*beta*c + 2*Q*beta*co - Q*beta*ct - Q*beta*w - 108*Q^2*delta^2*gamma^2 - 58*Q^2*delta^2*gamma^3 - 85*Q^2*delta^3*gamma^2 - 11*Q^2*delta^2*gamma^4 - 23*Q^2*delta^3*gamma^3 - 23*Q^2*delta^4*gamma^2 - 57*Q^2*delta^2*k^2 - 97*Q^2*delta^3*k^2 - 89*Q^2*delta^4*k^2 - 42*Q^2*delta^5*k^2 - 8*Q^2*delta^6*k^2 - 12*Q^2*gamma^2*k^2 - 8*Q^2*gamma^3*k^2 - 2*Q^2*gamma^4*k^2 + 14*beta^2*c^2*delta^2 + 37*beta^2*c^2*delta^3 + 47*beta^2*c^2*delta^4 + 28*beta^2*c^2*delta^5 + 6*beta^2*c^2*delta^6 - 25*beta^2*co^2*delta^2 - 39*beta^2*co^2*delta^3 - 32*beta^2*co^2*delta^4 - 13*beta^2*co^2*delta^5 - 2*beta^2*co^2*delta^6 + 7*beta^2*ct^2*delta^2 + 18*beta^2*ct^2*delta^3 + 22*beta^2*ct^2*delta^4 + 13*beta^2*ct^2*delta^5 + 3*beta^2*ct^2*delta^6 - 10*beta^2*co^2*gamma^2 - 10*beta^2*co^2*gamma^3 - 5*beta^2*co^2*gamma^4 - beta^2*co^2*gamma^5 + 4*beta^2*ct^2*gamma^2 + 6*beta^2*ct^2*gamma^3 + 4*beta^2*ct^2*gamma^4 + beta^2*ct^2*gamma^5 + 89*beta^2*delta^2*w^2 + 194*beta^2*delta^3*w^2 + 228*beta^2*delta^4*w^2 + 136*beta^2*delta^5*w^2 + 32*beta^2*delta^6*w^2 + 12*beta^2*gamma^2*w^2 + 8*beta^2*gamma^3*w^2 + 2*beta^2*gamma^4*w^2 - 34*Q^2*delta*gamma + 26*Q^2*delta*k + 14*Q^2*gamma*k - 2*beta^2*c*co + beta^2*co*ct - 2*beta^2*c*w + beta^2*co*w + beta^2*ct*w - 72*Q^2*delta*gamma^2*k^2 - 158*Q^2*delta^2*gamma*k^2 + 348*Q^2*delta^2*gamma^2*k - 38*Q^2*delta*gamma^3*k^2 + 176*Q^2*delta^2*gamma^3*k - 204*Q^2*delta^3*gamma*k^2 + 314*Q^2*delta^3*gamma^2*k - 7*Q^2*delta*gamma^4*k^2 + 31*Q^2*delta^2*gamma^4*k + 82*Q^2*delta^3*gamma^3*k - 126*Q^2*delta^4*gamma*k^2 + 100*Q^2*delta^4*gamma^2*k - 30*Q^2*delta^5*gamma*k^2 + 6*beta^2*c^2*delta*gamma^2 + 30*beta^2*c^2*delta^2*gamma + 2*beta^2*c^2*delta*gamma^3 + 46*beta^2*c^2*delta^3*gamma + 26*beta^2*c^2*delta^4*gamma + 4*beta^2*c^2*delta^5*gamma - 56*beta^2*co^2*delta*gamma^2 - 86*beta^2*co^2*delta^2*gamma - 44*beta^2*co^2*delta*gamma^3 - 101*beta^2*co^2*delta^3*gamma - 16*beta^2*co^2*delta*gamma^4 - 55*beta^2*co^2*delta^4*gamma - 2*beta^2*co^2*delta*gamma^5 - 11*beta^2*co^2*delta^5*gamma + 32*beta^2*ct^2*delta*gamma^2 + 48*beta^2*ct^2*delta^2*gamma + 34*beta^2*ct^2*delta*gamma^3 + 81*beta^2*ct^2*delta^3*gamma + 15*beta^2*ct^2*delta*gamma^4 + 61*beta^2*ct^2*delta^4*gamma + 2*beta^2*ct^2*delta*gamma^5 + 17*beta^2*ct^2*delta^5*gamma + 84*beta^2*delta*gamma^2*w^2 + 236*beta^2*delta^2*gamma*w^2 + 42*beta^2*delta*gamma^3*w^2 + 382*beta^2*delta^3*gamma*w^2 + 7*beta^2*delta*gamma^4*w^2 + 296*beta^2*delta^4*gamma*w^2 + 88*beta^2*delta^5*gamma*w^2 + 59*Q*beta*c*delta^2 + 108*Q*beta*c*delta^3 + 110*Q*beta*c*delta^4 + 59*Q*beta*c*delta^5 + 13*Q*beta*c*delta^6 + 50*Q*beta*co*delta^2 + 78*Q*beta*co*delta^3 + 64*Q*beta*co*delta^4 + 26*Q*beta*co*delta^5 + 4*Q*beta*co*delta^6 - 39*Q*beta*ct*delta^2 - 76*Q*beta*ct*delta^3 - 79*Q*beta*ct*delta^4 - 42*Q*beta*ct*delta^5 - 9*Q*beta*ct*delta^6 + 20*Q*beta*c*gamma^2 + 20*Q*beta*c*gamma^3 + 10*Q*beta*c*gamma^4 + 2*Q*beta*c*gamma^5 + 20*Q*beta*co*gamma^2 + 20*Q*beta*co*gamma^3 + 10*Q*beta*co*gamma^4 + 2*Q*beta*co*gamma^5 - 6*Q*beta*ct*gamma^2 - 4*Q*beta*ct*gamma^3 - Q*beta*ct*gamma^4 + 104*Q^2*delta*gamma*k - 17*beta^2*c*co*delta + beta^2*c*ct*delta + 10*beta^2*co*ct*delta - 48*Q*beta*delta^2*w - 106*Q*beta*delta^3*w - 125*Q*beta*delta^4*w - 75*Q*beta*delta^5*w - 18*Q*beta*delta^6*w - 10*beta^2*c*co*gamma + 2*beta^2*c*ct*gamma + 4*beta^2*co*ct*gamma - 6*Q*beta*gamma^2*w - 4*Q*beta*gamma^3*w - Q*beta*gamma^4*w - 20*beta^2*c*delta*w + 11*beta^2*co*delta*w + 11*beta^2*ct*delta*w - 8*beta^2*c*gamma*w + 4*beta^2*co*gamma*w + 4*beta^2*ct*gamma*w - 150*Q^2*delta^2*gamma^2*k^2 - 54*Q^2*delta^2*gamma^3*k^2 - 131*Q^2*delta^3*gamma^2*k^2 - 5*Q^2*delta^2*gamma^4*k^2 - 24*Q^2*delta^3*gamma^3*k^2 - 41*Q^2*delta^4*gamma^2*k^2 + 18*beta^2*c^2*delta^2*gamma^2 + 2*beta^2*c^2*delta^2*gamma^3 + beta^2*c^2*delta^3*gamma^2 - 8*beta^2*c^2*delta^3*gamma^3 - 11*beta^2*c^2*delta^4*gamma^2 - 108*beta^2*co^2*delta^2*gamma^2 - 58*beta^2*co^2*delta^2*gamma^3 - 85*beta^2*co^2*delta^3*gamma^2 - 11*beta^2*co^2*delta^2*gamma^4 - 23*beta^2*co^2*delta^3*gamma^3 - 23*beta^2*co^2*delta^4*gamma^2 + 90*beta^2*ct^2*delta^2*gamma^2 + 64*beta^2*ct^2*delta^2*gamma^3 + 98*beta^2*ct^2*delta^3*gamma^2 + 15*beta^2*ct^2*delta^2*gamma^4 + 35*beta^2*ct^2*delta^3*gamma^3 + 36*beta^2*ct^2*delta^4*gamma^2 + 210*beta^2*delta^2*gamma^2*w^2 + 68*beta^2*delta^2*gamma^3*w^2 + 222*beta^2*delta^3*gamma^2*w^2 + 5*beta^2*delta^2*gamma^4*w^2 + 34*beta^2*delta^3*gamma^3*w^2 + 84*beta^2*delta^4*gamma^2*w^2 - 58*Q^2*delta*gamma*k^2 + 160*Q^2*delta*gamma^2*k + 292*Q^2*delta^2*gamma*k + 116*Q^2*delta*gamma^3*k + 386*Q^2*delta^3*gamma*k + 38*Q^2*delta*gamma^4*k + 242*Q^2*delta^4*gamma*k + 4*Q^2*delta*gamma^5*k + 58*Q^2*delta^5*gamma*k - 59*beta^2*c*co*delta^2 - 108*beta^2*c*co*delta^3 - 110*beta^2*c*co*delta^4 - 59*beta^2*c*co*delta^5 - 13*beta^2*c*co*delta^6 + 5*beta^2*c*ct*delta^2 + 6*beta^2*c*ct*delta^3 - 2*beta^2*c*ct*delta^4 - 7*beta^2*c*ct*delta^5 - 3*beta^2*c*ct*delta^6 + 39*beta^2*co*ct*delta^2 + 76*beta^2*co*ct*delta^3 + 79*beta^2*co*ct*delta^4 + 42*beta^2*co*ct*delta^5 + 9*beta^2*co*ct*delta^6 - 20*beta^2*c*co*gamma^2 - 20*beta^2*c*co*gamma^3 - 10*beta^2*c*co*gamma^4 - 2*beta^2*c*co*gamma^5 + 8*beta^2*c*ct*gamma^2 + 12*beta^2*c*ct*gamma^3 + 8*beta^2*c*ct*gamma^4 + 2*beta^2*c*ct*gamma^5 + 6*beta^2*co*ct*gamma^2 + 4*beta^2*co*ct*gamma^3 + beta^2*co*ct*gamma^4 + 6*beta^2*c^2*delta*gamma - 34*beta^2*co^2*delta*gamma + 12*beta^2*ct^2*delta*gamma - 82*beta^2*c*delta^2*w - 176*beta^2*c*delta^3*w - 206*beta^2*c*delta^4*w - 122*beta^2*c*delta^5*w - 28*beta^2*c*delta^6*w + 48*beta^2*co*delta^2*w + 106*beta^2*co*delta^3*w + 125*beta^2*co*delta^4*w + 75*beta^2*co*delta^5*w + 18*beta^2*co*delta^6*w + 48*beta^2*ct*delta^2*w + 106*beta^2*ct*delta^3*w + 125*beta^2*ct*delta^4*w + 75*beta^2*ct*delta^5*w + 18*beta^2*ct*delta^6*w - 12*beta^2*c*gamma^2*w - 8*beta^2*c*gamma^3*w - 2*beta^2*c*gamma^4*w + 6*beta^2*co*gamma^2*w + 4*beta^2*co*gamma^3*w + beta^2*co*gamma^4*w + 6*beta^2*ct*gamma^2*w + 4*beta^2*ct*gamma^3*w + beta^2*ct*gamma^4*w + 70*beta^2*delta*gamma*w^2 + 17*Q*beta*c*delta + 16*Q*beta*co*delta - 10*Q*beta*ct*delta + 10*Q*beta*c*gamma + 10*Q*beta*co*gamma - 4*Q*beta*ct*gamma - 2*Q*beta*c*k - 3*Q*beta*co*k + Q*beta*ct*k - 11*Q*beta*delta*w - 4*Q*beta*gamma*w + 237*Q*beta*c*delta^2*gamma^2 + 122*Q*beta*c*delta^2*gamma^3 + 215*Q*beta*c*delta^3*gamma^2 + 22*Q*beta*c*delta^2*gamma^4 + 55*Q*beta*c*delta^3*gamma^3 + 73*Q*beta*c*delta^4*gamma^2 + 216*Q*beta*co*delta^2*gamma^2 + 116*Q*beta*co*delta^2*gamma^3 + 170*Q*beta*co*delta^3*gamma^2 + 22*Q*beta*co*delta^2*gamma^4 + 46*Q*beta*co*delta^3*gamma^3 + 46*Q*beta*co*delta^4*gamma^2 - 132*Q*beta*ct*delta^2*gamma^2 - 60*Q*beta*ct*delta^2*gamma^3 - 144*Q*beta*ct*delta^3*gamma^2 - 9*Q*beta*ct*delta^2*gamma^4 - 36*Q*beta*ct*delta^3*gamma^3 - 54*Q*beta*ct*delta^4*gamma^2 - 115*beta^2*c*co*delta*gamma^2 - 196*beta^2*c*co*delta^2*gamma - 89*beta^2*c*co*delta*gamma^3 - 268*beta^2*c*co*delta^3*gamma - 32*beta^2*c*co*delta*gamma^4 - 182*beta^2*c*co*delta^4*gamma - 4*beta^2*c*co*delta*gamma^5 - 49*beta^2*c*co*delta^5*gamma + 61*beta^2*c*ct*delta*gamma^2 + 72*beta^2*c*ct*delta^2*gamma + 67*beta^2*c*ct*delta*gamma^3 + 96*beta^2*c*ct*delta^3*gamma + 30*beta^2*c*ct*delta*gamma^4 + 50*beta^2*c*ct*delta^4*gamma + 4*beta^2*c*ct*delta*gamma^5 + 7*beta^2*c*ct*delta^5*gamma + 48*beta^2*co*ct*delta*gamma^2 + 120*beta^2*co*ct*delta^2*gamma + 28*beta^2*co*ct*delta*gamma^3 + 184*beta^2*co*ct*delta^3*gamma + 6*beta^2*co*ct*delta*gamma^4 + 132*beta^2*co*ct*delta^4*gamma + 36*beta^2*co*ct*delta^5*gamma - 153*Q*beta*delta^2*gamma^2*w - 66*Q*beta*delta^2*gamma^3*w - 189*Q*beta*delta^3*gamma^2*w - 9*Q*beta*delta^2*gamma^4*w - 45*Q*beta*delta^3*gamma^3*w - 81*Q*beta*delta^4*gamma^2*w - 66*beta^2*c*delta*gamma^2*w - 184*beta^2*c*delta^2*gamma*w - 26*beta^2*c*delta*gamma^3*w - 264*beta^2*c*delta^3*gamma*w - 2*beta^2*c*delta*gamma^4*w - 184*beta^2*c*delta^4*gamma*w - 50*beta^2*c*delta^5*gamma*w + 51*beta^2*co*delta*gamma^2*w + 144*beta^2*co*delta^2*gamma*w + 29*beta^2*co*delta*gamma^3*w + 250*beta^2*co*delta^3*gamma*w + 6*beta^2*co*delta*gamma^4*w + 204*beta^2*co*delta^4*gamma*w + 63*beta^2*co*delta^5*gamma*w + 51*beta^2*ct*delta*gamma^2*w + 144*beta^2*ct*delta^2*gamma*w + 29*beta^2*ct*delta*gamma^3*w + 250*beta^2*ct*delta^3*gamma*w + 6*beta^2*ct*delta*gamma^4*w + 204*beta^2*ct*delta^4*gamma*w + 63*beta^2*ct*delta^5*gamma*w + 71*Q*beta*c*delta*gamma + 68*Q*beta*co*delta*gamma - 36*Q*beta*ct*delta*gamma - 18*Q*beta*c*delta*k - 26*Q*beta*co*delta*k + 8*Q*beta*ct*delta*k - 12*Q*beta*c*gamma*k - 14*Q*beta*co*gamma*k + 2*Q*beta*ct*gamma*k - 39*Q*beta*delta*gamma*w - 237*beta^2*c*co*delta^2*gamma^2 - 122*beta^2*c*co*delta^2*gamma^3 - 215*beta^2*c*co*delta^3*gamma^2 - 22*beta^2*c*co*delta^2*gamma^4 - 55*beta^2*c*co*delta^3*gamma^3 - 73*beta^2*c*co*delta^4*gamma^2 + 159*beta^2*c*ct*delta^2*gamma^2 + 122*beta^2*c*ct*delta^2*gamma^3 + 151*beta^2*c*ct*delta^3*gamma^2 + 30*beta^2*c*ct*delta^2*gamma^4 + 61*beta^2*c*ct*delta^3*gamma^3 + 45*beta^2*c*ct*delta^4*gamma^2 + 132*beta^2*co*ct*delta^2*gamma^2 + 60*beta^2*co*ct*delta^2*gamma^3 + 144*beta^2*co*ct*delta^3*gamma^2 + 9*beta^2*co*ct*delta^2*gamma^4 + 36*beta^2*co*ct*delta^3*gamma^3 + 54*beta^2*co*ct*delta^4*gamma^2 - 114*beta^2*c*delta^2*gamma^2*w - 4*beta^2*c*delta^2*gamma^3*w - 66*beta^2*c*delta^3*gamma^2*w + 8*beta^2*c*delta^2*gamma^4*w + 22*beta^2*c*delta^3*gamma^3*w - 6*beta^2*c*delta^4*gamma^2*w + 153*beta^2*co*delta^2*gamma^2*w + 66*beta^2*co*delta^2*gamma^3*w + 189*beta^2*co*delta^3*gamma^2*w + 9*beta^2*co*delta^2*gamma^4*w + 45*beta^2*co*delta^3*gamma^3*w + 81*beta^2*co*delta^4*gamma^2*w + 153*beta^2*ct*delta^2*gamma^2*w + 66*beta^2*ct*delta^2*gamma^3*w + 189*beta^2*ct*delta^3*gamma^2*w + 9*beta^2*ct*delta^2*gamma^4*w + 45*beta^2*ct*delta^3*gamma^3*w + 81*beta^2*ct*delta^4*gamma^2*w + 115*Q*beta*c*delta*gamma^2 + 196*Q*beta*c*delta^2*gamma + 89*Q*beta*c*delta*gamma^3 + 268*Q*beta*c*delta^3*gamma + 32*Q*beta*c*delta*gamma^4 + 182*Q*beta*c*delta^4*gamma + 4*Q*beta*c*delta*gamma^5 + 49*Q*beta*c*delta^5*gamma + 112*Q*beta*co*delta*gamma^2 + 172*Q*beta*co*delta^2*gamma + 88*Q*beta*co*delta*gamma^3 + 202*Q*beta*co*delta^3*gamma + 32*Q*beta*co*delta*gamma^4 + 110*Q*beta*co*delta^4*gamma + 4*Q*beta*co*delta*gamma^5 + 22*Q*beta*co*delta^5*gamma - 48*Q*beta*ct*delta*gamma^2 - 120*Q*beta*ct*delta^2*gamma - 28*Q*beta*ct*delta*gamma^3 - 184*Q*beta*ct*delta^3*gamma - 6*Q*beta*ct*delta*gamma^4 - 132*Q*beta*ct*delta^4*gamma - 36*Q*beta*ct*delta^5*gamma - 64*Q*beta*c*delta^2*k - 114*Q*beta*c*delta^3*k - 108*Q*beta*c*delta^4*k - 52*Q*beta*c*delta^5*k - 10*Q*beta*c*delta^6*k - 89*Q*beta*co*delta^2*k - 154*Q*beta*co*delta^3*k - 143*Q*beta*co*delta^4*k - 68*Q*beta*co*delta^5*k - 13*Q*beta*co*delta^6*k + 25*Q*beta*ct*delta^2*k + 40*Q*beta*ct*delta^3*k + 35*Q*beta*ct*delta^4*k + 16*Q*beta*ct*delta^5*k + 3*Q*beta*ct*delta^6*k - 28*Q*beta*c*gamma^2*k - 32*Q*beta*c*gamma^3*k - 18*Q*beta*c*gamma^4*k - 4*Q*beta*c*gamma^5*k - 26*Q*beta*co*gamma^2*k - 24*Q*beta*co*gamma^3*k - 11*Q*beta*co*gamma^4*k - 2*Q*beta*co*gamma^5*k - 2*Q*beta*ct*gamma^2*k - 8*Q*beta*ct*gamma^3*k - 7*Q*beta*ct*gamma^4*k - 2*Q*beta*ct*gamma^5*k - 71*beta^2*c*co*delta*gamma + 21*beta^2*c*ct*delta*gamma + 36*beta^2*co*ct*delta*gamma - 51*Q*beta*delta*gamma^2*w - 144*Q*beta*delta^2*gamma*w - 29*Q*beta*delta*gamma^3*w - 250*Q*beta*delta^3*gamma*w - 6*Q*beta*delta*gamma^4*w - 204*Q*beta*delta^4*gamma*w - 63*Q*beta*delta^5*gamma*w - 62*beta^2*c*delta*gamma*w + 39*beta^2*co*delta*gamma*w + 39*beta^2*ct*delta*gamma*w - 176*Q*beta*c*delta*gamma^2*k - 268*Q*beta*c*delta^2*gamma*k - 156*Q*beta*c*delta*gamma^3*k - 364*Q*beta*c*delta^3*gamma*k - 62*Q*beta*c*delta*gamma^4*k - 232*Q*beta*c*delta^4*gamma*k - 8*Q*beta*c*delta*gamma^5*k - 56*Q*beta*c*delta^5*gamma*k - 160*Q*beta*co*delta*gamma^2*k - 292*Q*beta*co*delta^2*gamma*k - 116*Q*beta*co*delta*gamma^3*k - 386*Q*beta*co*delta^3*gamma*k - 38*Q*beta*co*delta*gamma^4*k - 242*Q*beta*co*delta^4*gamma*k - 4*Q*beta*co*delta*gamma^5*k - 58*Q*beta*co*delta^5*gamma*k - 16*Q*beta*ct*delta*gamma^2*k + 24*Q*beta*ct*delta^2*gamma*k - 40*Q*beta*ct*delta*gamma^3*k + 22*Q*beta*ct*delta^3*gamma*k - 24*Q*beta*ct*delta*gamma^4*k + 10*Q*beta*ct*delta^4*gamma*k - 4*Q*beta*ct*delta*gamma^5*k + 2*Q*beta*ct*delta^5*gamma*k - 396*Q*beta*c*delta^2*gamma^2*k - 244*Q*beta*c*delta^2*gamma^3*k - 366*Q*beta*c*delta^3*gamma^2*k - 52*Q*beta*c*delta^2*gamma^4*k - 116*Q*beta*c*delta^3*gamma^3*k - 118*Q*beta*c*delta^4*gamma^2*k - 348*Q*beta*co*delta^2*gamma^2*k - 176*Q*beta*co*delta^2*gamma^3*k - 314*Q*beta*co*delta^3*gamma^2*k - 31*Q*beta*co*delta^2*gamma^4*k - 82*Q*beta*co*delta^3*gamma^3*k - 100*Q*beta*co*delta^4*gamma^2*k - 48*Q*beta*ct*delta^2*gamma^2*k - 68*Q*beta*ct*delta^2*gamma^3*k - 52*Q*beta*ct*delta^3*gamma^2*k - 21*Q*beta*ct*delta^2*gamma^4*k - 34*Q*beta*ct*delta^3*gamma^3*k - 18*Q*beta*ct*delta^4*gamma^2*k - 92*Q*beta*c*delta*gamma*k - 104*Q*beta*co*delta*gamma*k + 12*Q*beta*ct*delta*gamma*k)/(4*beta*(delta + gamma + 1)^2*(11*delta + 4*gamma + 16*delta*gamma + 5*delta*gamma^2 + 14*delta^2*gamma + 18*delta^2 + 8*delta^3 + 2*gamma^2 + 2));
Utn=-(3*Q^2*delta^6*k - 8*Q^2*delta^6*k^2 + 3*Q^2*delta^6 - 30*Q^2*delta^5*gamma*k^2 + 2*Q^2*delta^5*gamma*k + 17*Q^2*delta^5*gamma - 42*Q^2*delta^5*k^2 + 16*Q^2*delta^5*k + 13*Q^2*delta^5 - 41*Q^2*delta^4*gamma^2*k^2 - 18*Q^2*delta^4*gamma^2*k + 36*Q^2*delta^4*gamma^2 - 126*Q^2*delta^4*gamma*k^2 + 10*Q^2*delta^4*gamma*k + 61*Q^2*delta^4*gamma - 89*Q^2*delta^4*k^2 + 35*Q^2*delta^4*k + 22*Q^2*delta^4 - 24*Q^2*delta^3*gamma^3*k^2 - 34*Q^2*delta^3*gamma^3*k + 35*Q^2*delta^3*gamma^3 - 131*Q^2*delta^3*gamma^2*k^2 - 52*Q^2*delta^3*gamma^2*k + 98*Q^2*delta^3*gamma^2 - 204*Q^2*delta^3*gamma*k^2 + 22*Q^2*delta^3*gamma*k + 81*Q^2*delta^3*gamma - 97*Q^2*delta^3*k^2 + 40*Q^2*delta^3*k + 18*Q^2*delta^3 - 5*Q^2*delta^2*gamma^4*k^2 - 21*Q^2*delta^2*gamma^4*k + 15*Q^2*delta^2*gamma^4 - 54*Q^2*delta^2*gamma^3*k^2 - 68*Q^2*delta^2*gamma^3*k + 64*Q^2*delta^2*gamma^3 - 150*Q^2*delta^2*gamma^2*k^2 - 48*Q^2*delta^2*gamma^2*k + 90*Q^2*delta^2*gamma^2 - 158*Q^2*delta^2*gamma*k^2 + 24*Q^2*delta^2*gamma*k + 48*Q^2*delta^2*gamma - 57*Q^2*delta^2*k^2 + 25*Q^2*delta^2*k + 7*Q^2*delta^2 - 4*Q^2*delta*gamma^5*k + 2*Q^2*delta*gamma^5 - 7*Q^2*delta*gamma^4*k^2 - 24*Q^2*delta*gamma^4*k + 15*Q^2*delta*gamma^4 - 38*Q^2*delta*gamma^3*k^2 - 40*Q^2*delta*gamma^3*k + 34*Q^2*delta*gamma^3 - 72*Q^2*delta*gamma^2*k^2 - 16*Q^2*delta*gamma^2*k + 32*Q^2*delta*gamma^2 - 58*Q^2*delta*gamma*k^2 + 12*Q^2*delta*gamma*k + 12*Q^2*delta*gamma - 17*Q^2*delta*k^2 + 8*Q^2*delta*k + Q^2*delta - 2*Q^2*gamma^5*k + Q^2*gamma^5 - 2*Q^2*gamma^4*k^2 - 7*Q^2*gamma^4*k + 4*Q^2*gamma^4 - 8*Q^2*gamma^3*k^2 - 8*Q^2*gamma^3*k + 6*Q^2*gamma^3 - 12*Q^2*gamma^2*k^2 - 2*Q^2*gamma^2*k + 4*Q^2*gamma^2 - 8*Q^2*gamma*k^2 + 2*Q^2*gamma*k + Q^2*gamma - 2*Q^2*k^2 + Q^2*k + 10*Q*beta*c*delta^6*k + 3*Q*beta*c*delta^6 + 56*Q*beta*c*delta^5*gamma*k - 7*Q*beta*c*delta^5*gamma + 52*Q*beta*c*delta^5*k + 7*Q*beta*c*delta^5 + 118*Q*beta*c*delta^4*gamma^2*k - 45*Q*beta*c*delta^4*gamma^2 + 232*Q*beta*c*delta^4*gamma*k - 50*Q*beta*c*delta^4*gamma + 108*Q*beta*c*delta^4*k + 2*Q*beta*c*delta^4 + 116*Q*beta*c*delta^3*gamma^3*k - 61*Q*beta*c*delta^3*gamma^3 + 366*Q*beta*c*delta^3*gamma^2*k - 151*Q*beta*c*delta^3*gamma^2 + 364*Q*beta*c*delta^3*gamma*k - 96*Q*beta*c*delta^3*gamma + 114*Q*beta*c*delta^3*k - 6*Q*beta*c*delta^3 + 52*Q*beta*c*delta^2*gamma^4*k - 30*Q*beta*c*delta^2*gamma^4 + 244*Q*beta*c*delta^2*gamma^3*k - 122*Q*beta*c*delta^2*gamma^3 + 396*Q*beta*c*delta^2*gamma^2*k - 159*Q*beta*c*delta^2*gamma^2 + 268*Q*beta*c*delta^2*gamma*k - 72*Q*beta*c*delta^2*gamma + 64*Q*beta*c*delta^2*k - 5*Q*beta*c*delta^2 + 8*Q*beta*c*delta*gamma^5*k - 4*Q*beta*c*delta*gamma^5 + 62*Q*beta*c*delta*gamma^4*k - 30*Q*beta*c*delta*gamma^4 + 156*Q*beta*c*delta*gamma^3*k - 67*Q*beta*c*delta*gamma^3 + 176*Q*beta*c*delta*gamma^2*k - 61*Q*beta*c*delta*gamma^2 + 92*Q*beta*c*delta*gamma*k - 21*Q*beta*c*delta*gamma + 18*Q*beta*c*delta*k - Q*beta*c*delta + 4*Q*beta*c*gamma^5*k - 2*Q*beta*c*gamma^5 + 18*Q*beta*c*gamma^4*k - 8*Q*beta*c*gamma^4 + 32*Q*beta*c*gamma^3*k - 12*Q*beta*c*gamma^3 + 28*Q*beta*c*gamma^2*k - 8*Q*beta*c*gamma^2 + 12*Q*beta*c*gamma*k - 2*Q*beta*c*gamma + 2*Q*beta*c*k - 3*Q*beta*co*delta^6*k - 6*Q*beta*co*delta^6 - 2*Q*beta*co*delta^5*gamma*k - 34*Q*beta*co*delta^5*gamma - 16*Q*beta*co*delta^5*k - 26*Q*beta*co*delta^5 + 18*Q*beta*co*delta^4*gamma^2*k - 72*Q*beta*co*delta^4*gamma^2 - 10*Q*beta*co*delta^4*gamma*k - 122*Q*beta*co*delta^4*gamma - 35*Q*beta*co*delta^4*k - 44*Q*beta*co*delta^4 + 34*Q*beta*co*delta^3*gamma^3*k - 70*Q*beta*co*delta^3*gamma^3 + 52*Q*beta*co*delta^3*gamma^2*k - 196*Q*beta*co*delta^3*gamma^2 - 22*Q*beta*co*delta^3*gamma*k - 162*Q*beta*co*delta^3*gamma - 40*Q*beta*co*delta^3*k - 36*Q*beta*co*delta^3 + 21*Q*beta*co*delta^2*gamma^4*k - 30*Q*beta*co*delta^2*gamma^4 + 68*Q*beta*co*delta^2*gamma^3*k - 128*Q*beta*co*delta^2*gamma^3 + 48*Q*beta*co*delta^2*gamma^2*k - 180*Q*beta*co*delta^2*gamma^2 - 24*Q*beta*co*delta^2*gamma*k - 96*Q*beta*co*delta^2*gamma - 25*Q*beta*co*delta^2*k - 14*Q*beta*co*delta^2 + 4*Q*beta*co*delta*gamma^5*k - 4*Q*beta*co*delta*gamma^5 + 24*Q*beta*co*delta*gamma^4*k - 30*Q*beta*co*delta*gamma^4 + 40*Q*beta*co*delta*gamma^3*k - 68*Q*beta*co*delta*gamma^3 + 16*Q*beta*co*delta*gamma^2*k - 64*Q*beta*co*delta*gamma^2 - 12*Q*beta*co*delta*gamma*k - 24*Q*beta*co*delta*gamma - 8*Q*beta*co*delta*k - 2*Q*beta*co*delta + 2*Q*beta*co*gamma^5*k - 2*Q*beta*co*gamma^5 + 7*Q*beta*co*gamma^4*k - 8*Q*beta*co*gamma^4 + 8*Q*beta*co*gamma^3*k - 12*Q*beta*co*gamma^3 + 2*Q*beta*co*gamma^2*k - 8*Q*beta*co*gamma^2 - 2*Q*beta*co*gamma*k - 2*Q*beta*co*gamma - Q*beta*co*k + 13*Q*beta*ct*delta^6*k - 9*Q*beta*ct*delta^6 + 58*Q*beta*ct*delta^5*gamma*k - 36*Q*beta*ct*delta^5*gamma + 68*Q*beta*ct*delta^5*k - 42*Q*beta*ct*delta^5 + 100*Q*beta*ct*delta^4*gamma^2*k - 54*Q*beta*ct*delta^4*gamma^2 + 242*Q*beta*ct*delta^4*gamma*k - 132*Q*beta*ct*delta^4*gamma + 143*Q*beta*ct*delta^4*k - 79*Q*beta*ct*delta^4 + 82*Q*beta*ct*delta^3*gamma^3*k - 36*Q*beta*ct*delta^3*gamma^3 + 314*Q*beta*ct*delta^3*gamma^2*k - 144*Q*beta*ct*delta^3*gamma^2 + 386*Q*beta*ct*delta^3*gamma*k - 184*Q*beta*ct*delta^3*gamma + 154*Q*beta*ct*delta^3*k - 76*Q*beta*ct*delta^3 + 31*Q*beta*ct*delta^2*gamma^4*k - 9*Q*beta*ct*delta^2*gamma^4 + 176*Q*beta*ct*delta^2*gamma^3*k - 60*Q*beta*ct*delta^2*gamma^3 + 348*Q*beta*ct*delta^2*gamma^2*k - 132*Q*beta*ct*delta^2*gamma^2 + 292*Q*beta*ct*delta^2*gamma*k - 120*Q*beta*ct*delta^2*gamma + 89*Q*beta*ct*delta^2*k - 39*Q*beta*ct*delta^2 + 4*Q*beta*ct*delta*gamma^5*k + 38*Q*beta*ct*delta*gamma^4*k - 6*Q*beta*ct*delta*gamma^4 + 116*Q*beta*ct*delta*gamma^3*k - 28*Q*beta*ct*delta*gamma^3 + 160*Q*beta*ct*delta*gamma^2*k - 48*Q*beta*ct*delta*gamma^2 + 104*Q*beta*ct*delta*gamma*k - 36*Q*beta*ct*delta*gamma + 26*Q*beta*ct*delta*k - 10*Q*beta*ct*delta + 2*Q*beta*ct*gamma^5*k + 11*Q*beta*ct*gamma^4*k - Q*beta*ct*gamma^4 + 24*Q*beta*ct*gamma^3*k - 4*Q*beta*ct*gamma^3 + 26*Q*beta*ct*gamma^2*k - 6*Q*beta*ct*gamma^2 + 14*Q*beta*ct*gamma*k - 4*Q*beta*ct*gamma + 3*Q*beta*ct*k - Q*beta*ct - 18*Q*beta*delta^6*w - 63*Q*beta*delta^5*gamma*w - 75*Q*beta*delta^5*w - 81*Q*beta*delta^4*gamma^2*w - 204*Q*beta*delta^4*gamma*w - 125*Q*beta*delta^4*w - 45*Q*beta*delta^3*gamma^3*w - 189*Q*beta*delta^3*gamma^2*w - 250*Q*beta*delta^3*gamma*w - 106*Q*beta*delta^3*w - 9*Q*beta*delta^2*gamma^4*w - 66*Q*beta*delta^2*gamma^3*w - 153*Q*beta*delta^2*gamma^2*w - 144*Q*beta*delta^2*gamma*w - 48*Q*beta*delta^2*w - 6*Q*beta*delta*gamma^4*w - 29*Q*beta*delta*gamma^3*w - 51*Q*beta*delta*gamma^2*w - 39*Q*beta*delta*gamma*w - 11*Q*beta*delta*w - Q*beta*gamma^4*w - 4*Q*beta*gamma^3*w - 6*Q*beta*gamma^2*w - 4*Q*beta*gamma*w - Q*beta*w + 6*beta^2*c^2*delta^6 + 4*beta^2*c^2*delta^5*gamma + 28*beta^2*c^2*delta^5 - 11*beta^2*c^2*delta^4*gamma^2 + 26*beta^2*c^2*delta^4*gamma + 47*beta^2*c^2*delta^4 - 8*beta^2*c^2*delta^3*gamma^3 + beta^2*c^2*delta^3*gamma^2 + 46*beta^2*c^2*delta^3*gamma + 37*beta^2*c^2*delta^3 + 2*beta^2*c^2*delta^2*gamma^3 + 18*beta^2*c^2*delta^2*gamma^2 + 30*beta^2*c^2*delta^2*gamma + 14*beta^2*c^2*delta^2 + 2*beta^2*c^2*delta*gamma^3 + 6*beta^2*c^2*delta*gamma^2 + 6*beta^2*c^2*delta*gamma + 2*beta^2*c^2*delta - 3*beta^2*c*co*delta^6 + 7*beta^2*c*co*delta^5*gamma - 7*beta^2*c*co*delta^5 + 45*beta^2*c*co*delta^4*gamma^2 + 50*beta^2*c*co*delta^4*gamma - 2*beta^2*c*co*delta^4 + 61*beta^2*c*co*delta^3*gamma^3 + 151*beta^2*c*co*delta^3*gamma^2 + 96*beta^2*c*co*delta^3*gamma + 6*beta^2*c*co*delta^3 + 30*beta^2*c*co*delta^2*gamma^4 + 122*beta^2*c*co*delta^2*gamma^3 + 159*beta^2*c*co*delta^2*gamma^2 + 72*beta^2*c*co*delta^2*gamma + 5*beta^2*c*co*delta^2 + 4*beta^2*c*co*delta*gamma^5 + 30*beta^2*c*co*delta*gamma^4 + 67*beta^2*c*co*delta*gamma^3 + 61*beta^2*c*co*delta*gamma^2 + 21*beta^2*c*co*delta*gamma + beta^2*c*co*delta + 2*beta^2*c*co*gamma^5 + 8*beta^2*c*co*gamma^4 + 12*beta^2*c*co*gamma^3 + 8*beta^2*c*co*gamma^2 + 2*beta^2*c*co*gamma - 13*beta^2*c*ct*delta^6 - 49*beta^2*c*ct*delta^5*gamma - 59*beta^2*c*ct*delta^5 - 73*beta^2*c*ct*delta^4*gamma^2 - 182*beta^2*c*ct*delta^4*gamma - 110*beta^2*c*ct*delta^4 - 55*beta^2*c*ct*delta^3*gamma^3 - 215*beta^2*c*ct*delta^3*gamma^2 - 268*beta^2*c*ct*delta^3*gamma - 108*beta^2*c*ct*delta^3 - 22*beta^2*c*ct*delta^2*gamma^4 - 122*beta^2*c*ct*delta^2*gamma^3 - 237*beta^2*c*ct*delta^2*gamma^2 - 196*beta^2*c*ct*delta^2*gamma - 59*beta^2*c*ct*delta^2 - 4*beta^2*c*ct*delta*gamma^5 - 32*beta^2*c*ct*delta*gamma^4 - 89*beta^2*c*ct*delta*gamma^3 - 115*beta^2*c*ct*delta*gamma^2 - 71*beta^2*c*ct*delta*gamma - 17*beta^2*c*ct*delta - 2*beta^2*c*ct*gamma^5 - 10*beta^2*c*ct*gamma^4 - 20*beta^2*c*ct*gamma^3 - 20*beta^2*c*ct*gamma^2 - 10*beta^2*c*ct*gamma - 2*beta^2*c*ct - 28*beta^2*c*delta^6*w - 50*beta^2*c*delta^5*gamma*w - 122*beta^2*c*delta^5*w - 6*beta^2*c*delta^4*gamma^2*w - 184*beta^2*c*delta^4*gamma*w - 206*beta^2*c*delta^4*w + 22*beta^2*c*delta^3*gamma^3*w - 66*beta^2*c*delta^3*gamma^2*w - 264*beta^2*c*delta^3*gamma*w - 176*beta^2*c*delta^3*w + 8*beta^2*c*delta^2*gamma^4*w - 4*beta^2*c*delta^2*gamma^3*w - 114*beta^2*c*delta^2*gamma^2*w - 184*beta^2*c*delta^2*gamma*w - 82*beta^2*c*delta^2*w - 2*beta^2*c*delta*gamma^4*w - 26*beta^2*c*delta*gamma^3*w - 66*beta^2*c*delta*gamma^2*w - 62*beta^2*c*delta*gamma*w - 20*beta^2*c*delta*w - 2*beta^2*c*gamma^4*w - 8*beta^2*c*gamma^3*w - 12*beta^2*c*gamma^2*w - 8*beta^2*c*gamma*w - 2*beta^2*c*w + 3*beta^2*co^2*delta^6 + 17*beta^2*co^2*delta^5*gamma + 13*beta^2*co^2*delta^5 + 36*beta^2*co^2*delta^4*gamma^2 + 61*beta^2*co^2*delta^4*gamma + 22*beta^2*co^2*delta^4 + 35*beta^2*co^2*delta^3*gamma^3 + 98*beta^2*co^2*delta^3*gamma^2 + 81*beta^2*co^2*delta^3*gamma + 18*beta^2*co^2*delta^3 + 15*beta^2*co^2*delta^2*gamma^4 + 64*beta^2*co^2*delta^2*gamma^3 + 90*beta^2*co^2*delta^2*gamma^2 + 48*beta^2*co^2*delta^2*gamma + 7*beta^2*co^2*delta^2 + 2*beta^2*co^2*delta*gamma^5 + 15*beta^2*co^2*delta*gamma^4 + 34*beta^2*co^2*delta*gamma^3 + 32*beta^2*co^2*delta*gamma^2 + 12*beta^2*co^2*delta*gamma + beta^2*co^2*delta + beta^2*co^2*gamma^5 + 4*beta^2*co^2*gamma^4 + 6*beta^2*co^2*gamma^3 + 4*beta^2*co^2*gamma^2 + beta^2*co^2*gamma + 9*beta^2*co*ct*delta^6 + 36*beta^2*co*ct*delta^5*gamma + 42*beta^2*co*ct*delta^5 + 54*beta^2*co*ct*delta^4*gamma^2 + 132*beta^2*co*ct*delta^4*gamma + 79*beta^2*co*ct*delta^4 + 36*beta^2*co*ct*delta^3*gamma^3 + 144*beta^2*co*ct*delta^3*gamma^2 + 184*beta^2*co*ct*delta^3*gamma + 76*beta^2*co*ct*delta^3 + 9*beta^2*co*ct*delta^2*gamma^4 + 60*beta^2*co*ct*delta^2*gamma^3 + 132*beta^2*co*ct*delta^2*gamma^2 + 120*beta^2*co*ct*delta^2*gamma + 39*beta^2*co*ct*delta^2 + 6*beta^2*co*ct*delta*gamma^4 + 28*beta^2*co*ct*delta*gamma^3 + 48*beta^2*co*ct*delta*gamma^2 + 36*beta^2*co*ct*delta*gamma + 10*beta^2*co*ct*delta + beta^2*co*ct*gamma^4 + 4*beta^2*co*ct*gamma^3 + 6*beta^2*co*ct*gamma^2 + 4*beta^2*co*ct*gamma + beta^2*co*ct + 18*beta^2*co*delta^6*w + 63*beta^2*co*delta^5*gamma*w + 75*beta^2*co*delta^5*w + 81*beta^2*co*delta^4*gamma^2*w + 204*beta^2*co*delta^4*gamma*w + 125*beta^2*co*delta^4*w + 45*beta^2*co*delta^3*gamma^3*w + 189*beta^2*co*delta^3*gamma^2*w + 250*beta^2*co*delta^3*gamma*w + 106*beta^2*co*delta^3*w + 9*beta^2*co*delta^2*gamma^4*w + 66*beta^2*co*delta^2*gamma^3*w + 153*beta^2*co*delta^2*gamma^2*w + 144*beta^2*co*delta^2*gamma*w + 48*beta^2*co*delta^2*w + 6*beta^2*co*delta*gamma^4*w + 29*beta^2*co*delta*gamma^3*w + 51*beta^2*co*delta*gamma^2*w + 39*beta^2*co*delta*gamma*w + 11*beta^2*co*delta*w + beta^2*co*gamma^4*w + 4*beta^2*co*gamma^3*w + 6*beta^2*co*gamma^2*w + 4*beta^2*co*gamma*w + beta^2*co*w - 2*beta^2*ct^2*delta^6 - 11*beta^2*ct^2*delta^5*gamma - 13*beta^2*ct^2*delta^5 - 23*beta^2*ct^2*delta^4*gamma^2 - 55*beta^2*ct^2*delta^4*gamma - 32*beta^2*ct^2*delta^4 - 23*beta^2*ct^2*delta^3*gamma^3 - 85*beta^2*ct^2*delta^3*gamma^2 - 101*beta^2*ct^2*delta^3*gamma - 39*beta^2*ct^2*delta^3 - 11*beta^2*ct^2*delta^2*gamma^4 - 58*beta^2*ct^2*delta^2*gamma^3 - 108*beta^2*ct^2*delta^2*gamma^2 - 86*beta^2*ct^2*delta^2*gamma - 25*beta^2*ct^2*delta^2 - 2*beta^2*ct^2*delta*gamma^5 - 16*beta^2*ct^2*delta*gamma^4 - 44*beta^2*ct^2*delta*gamma^3 - 56*beta^2*ct^2*delta*gamma^2 - 34*beta^2*ct^2*delta*gamma - 8*beta^2*ct^2*delta - beta^2*ct^2*gamma^5 - 5*beta^2*ct^2*gamma^4 - 10*beta^2*ct^2*gamma^3 - 10*beta^2*ct^2*gamma^2 - 5*beta^2*ct^2*gamma - beta^2*ct^2 + 18*beta^2*ct*delta^6*w + 63*beta^2*ct*delta^5*gamma*w + 75*beta^2*ct*delta^5*w + 81*beta^2*ct*delta^4*gamma^2*w + 204*beta^2*ct*delta^4*gamma*w + 125*beta^2*ct*delta^4*w + 45*beta^2*ct*delta^3*gamma^3*w + 189*beta^2*ct*delta^3*gamma^2*w + 250*beta^2*ct*delta^3*gamma*w + 106*beta^2*ct*delta^3*w + 9*beta^2*ct*delta^2*gamma^4*w + 66*beta^2*ct*delta^2*gamma^3*w + 153*beta^2*ct*delta^2*gamma^2*w + 144*beta^2*ct*delta^2*gamma*w + 48*beta^2*ct*delta^2*w + 6*beta^2*ct*delta*gamma^4*w + 29*beta^2*ct*delta*gamma^3*w + 51*beta^2*ct*delta*gamma^2*w + 39*beta^2*ct*delta*gamma*w + 11*beta^2*ct*delta*w + beta^2*ct*gamma^4*w + 4*beta^2*ct*gamma^3*w + 6*beta^2*ct*gamma^2*w + 4*beta^2*ct*gamma*w + beta^2*ct*w + 32*beta^2*delta^6*w^2 + 88*beta^2*delta^5*gamma*w^2 + 136*beta^2*delta^5*w^2 + 84*beta^2*delta^4*gamma^2*w^2 + 296*beta^2*delta^4*gamma*w^2 + 228*beta^2*delta^4*w^2 + 34*beta^2*delta^3*gamma^3*w^2 + 222*beta^2*delta^3*gamma^2*w^2 + 382*beta^2*delta^3*gamma*w^2 + 194*beta^2*delta^3*w^2 + 5*beta^2*delta^2*gamma^4*w^2 + 68*beta^2*delta^2*gamma^3*w^2 + 210*beta^2*delta^2*gamma^2*w^2 + 236*beta^2*delta^2*gamma*w^2 + 89*beta^2*delta^2*w^2 + 7*beta^2*delta*gamma^4*w^2 + 42*beta^2*delta*gamma^3*w^2 + 84*beta^2*delta*gamma^2*w^2 + 70*beta^2*delta*gamma*w^2 + 21*beta^2*delta*w^2 + 2*beta^2*gamma^4*w^2 + 8*beta^2*gamma^3*w^2 + 12*beta^2*gamma^2*w^2 + 8*beta^2*gamma*w^2 + 2*beta^2*w^2)/(4*beta*(delta + gamma + 1)^2*(8*delta^3 + 14*delta^2*gamma + 18*delta^2 + 5*delta*gamma^2 + 16*delta*gamma + 11*delta + 2*gamma^2 + 4*gamma + 2));
Un=(32*Q^2*delta^5*k^2 - 36*Q^2*delta^5*gamma - 9*Q^2*delta^6 - 32*Q^2*delta^5*k - 26*Q^2*delta^5 - 54*Q^2*delta^4*gamma^2 + 120*Q^2*delta^4*gamma*k^2 - 120*Q^2*delta^4*gamma*k - 72*Q^2*delta^4*gamma + 136*Q^2*delta^4*k^2 - 136*Q^2*delta^4*k - 11*Q^2*delta^4 - 36*Q^2*delta^3*gamma^3 + 164*Q^2*delta^3*gamma^2*k^2 - 164*Q^2*delta^3*gamma^2*k - 62*Q^2*delta^3*gamma^2 + 384*Q^2*delta^3*gamma*k^2 - 384*Q^2*delta^3*gamma*k + 8*Q^2*delta^3*gamma + 220*Q^2*delta^3*k^2 - 220*Q^2*delta^3*k + 34*Q^2*delta^3 - 9*Q^2*delta^2*gamma^4 + 96*Q^2*delta^2*gamma^3*k^2 - 96*Q^2*delta^2*gamma^3*k - 12*Q^2*delta^2*gamma^3 + 360*Q^2*delta^2*gamma^2*k^2 - 360*Q^2*delta^2*gamma^2*k + 48*Q^2*delta^2*gamma^2 + 432*Q^2*delta^2*gamma*k^2 - 432*Q^2*delta^2*gamma*k + 96*Q^2*delta^2*gamma + 168*Q^2*delta^2*k^2 - 168*Q^2*delta^2*k + 45*Q^2*delta^2 + 20*Q^2*delta*gamma^4*k^2 - 20*Q^2*delta*gamma^4*k + 4*Q^2*delta*gamma^4 + 120*Q^2*delta*gamma^3*k^2 - 120*Q^2*delta*gamma^3*k + 32*Q^2*delta*gamma^3 + 240*Q^2*delta*gamma^2*k^2 - 240*Q^2*delta*gamma^2*k + 72*Q^2*delta*gamma^2 + 200*Q^2*delta*gamma*k^2 - 200*Q^2*delta*gamma*k + 64*Q^2*delta*gamma + 60*Q^2*delta*k^2 - 60*Q^2*delta*k + 20*Q^2*delta + 8*Q^2*gamma^4*k^2 - 8*Q^2*gamma^4*k + 3*Q^2*gamma^4 + 32*Q^2*gamma^3*k^2 - 32*Q^2*gamma^3*k + 12*Q^2*gamma^3 + 48*Q^2*gamma^2*k^2 - 48*Q^2*gamma^2*k + 18*Q^2*gamma^2 + 32*Q^2*gamma*k^2 - 32*Q^2*gamma*k + 12*Q^2*gamma + 8*Q^2*k^2 - 8*Q^2*k + 3*Q^2 - 36*Q*beta*c*delta^6 - 108*Q*beta*c*delta^5*gamma - 164*Q*beta*c*delta^5 - 108*Q*beta*c*delta^4*gamma^2 - 408*Q*beta*c*delta^4*gamma - 320*Q*beta*c*delta^4 - 36*Q*beta*c*delta^3*gamma^3 - 344*Q*beta*c*delta^3*gamma^2 - 648*Q*beta*c*delta^3*gamma - 340*Q*beta*c*delta^3 - 120*Q*beta*c*delta^2*gamma^3 - 444*Q*beta*c*delta^2*gamma^2 - 528*Q*beta*c*delta^2*gamma - 204*Q*beta*c*delta^2 - 20*Q*beta*c*delta*gamma^4 - 124*Q*beta*c*delta*gamma^3 - 252*Q*beta*c*delta*gamma^2 - 212*Q*beta*c*delta*gamma - 64*Q*beta*c*delta - 8*Q*beta*c*gamma^4 - 32*Q*beta*c*gamma^3 - 48*Q*beta*c*gamma^2 - 32*Q*beta*c*gamma - 8*Q*beta*c + 18*Q*beta*co*delta^6 + 72*Q*beta*co*delta^5*gamma + 32*Q*beta*co*delta^5*k + 52*Q*beta*co*delta^5 + 108*Q*beta*co*delta^4*gamma^2 + 120*Q*beta*co*delta^4*gamma*k + 144*Q*beta*co*delta^4*gamma + 136*Q*beta*co*delta^4*k + 22*Q*beta*co*delta^4 + 72*Q*beta*co*delta^3*gamma^3 + 164*Q*beta*co*delta^3*gamma^2*k + 124*Q*beta*co*delta^3*gamma^2 + 384*Q*beta*co*delta^3*gamma*k - 16*Q*beta*co*delta^3*gamma + 220*Q*beta*co*delta^3*k - 68*Q*beta*co*delta^3 + 18*Q*beta*co*delta^2*gamma^4 + 96*Q*beta*co*delta^2*gamma^3*k + 24*Q*beta*co*delta^2*gamma^3 + 360*Q*beta*co*delta^2*gamma^2*k - 96*Q*beta*co*delta^2*gamma^2 + 432*Q*beta*co*delta^2*gamma*k - 192*Q*beta*co*delta^2*gamma + 168*Q*beta*co*delta^2*k - 90*Q*beta*co*delta^2 + 20*Q*beta*co*delta*gamma^4*k - 8*Q*beta*co*delta*gamma^4 + 120*Q*beta*co*delta*gamma^3*k - 64*Q*beta*co*delta*gamma^3 + 240*Q*beta*co*delta*gamma^2*k - 144*Q*beta*co*delta*gamma^2 + 200*Q*beta*co*delta*gamma*k - 128*Q*beta*co*delta*gamma + 60*Q*beta*co*delta*k - 40*Q*beta*co*delta + 8*Q*beta*co*gamma^4*k - 6*Q*beta*co*gamma^4 + 32*Q*beta*co*gamma^3*k - 24*Q*beta*co*gamma^3 + 48*Q*beta*co*gamma^2*k - 36*Q*beta*co*gamma^2 + 32*Q*beta*co*gamma*k - 24*Q*beta*co*gamma + 8*Q*beta*co*k - 6*Q*beta*co + 18*Q*beta*ct*delta^6 + 72*Q*beta*ct*delta^5*gamma - 32*Q*beta*ct*delta^5*k + 84*Q*beta*ct*delta^5 + 108*Q*beta*ct*delta^4*gamma^2 - 120*Q*beta*ct*delta^4*gamma*k + 264*Q*beta*ct*delta^4*gamma - 136*Q*beta*ct*delta^4*k + 158*Q*beta*ct*delta^4 + 72*Q*beta*ct*delta^3*gamma^3 - 164*Q*beta*ct*delta^3*gamma^2*k + 288*Q*beta*ct*delta^3*gamma^2 - 384*Q*beta*ct*delta^3*gamma*k + 368*Q*beta*ct*delta^3*gamma - 220*Q*beta*ct*delta^3*k + 152*Q*beta*ct*delta^3 + 18*Q*beta*ct*delta^2*gamma^4 - 96*Q*beta*ct*delta^2*gamma^3*k + 120*Q*beta*ct*delta^2*gamma^3 - 360*Q*beta*ct*delta^2*gamma^2*k + 264*Q*beta*ct*delta^2*gamma^2 - 432*Q*beta*ct*delta^2*gamma*k + 240*Q*beta*ct*delta^2*gamma - 168*Q*beta*ct*delta^2*k + 78*Q*beta*ct*delta^2 - 20*Q*beta*ct*delta*gamma^4*k + 12*Q*beta*ct*delta*gamma^4 - 120*Q*beta*ct*delta*gamma^3*k + 56*Q*beta*ct*delta*gamma^3 - 240*Q*beta*ct*delta*gamma^2*k + 96*Q*beta*ct*delta*gamma^2 - 200*Q*beta*ct*delta*gamma*k + 72*Q*beta*ct*delta*gamma - 60*Q*beta*ct*delta*k + 20*Q*beta*ct*delta - 8*Q*beta*ct*gamma^4*k + 2*Q*beta*ct*gamma^4 - 32*Q*beta*ct*gamma^3*k + 8*Q*beta*ct*gamma^3 - 48*Q*beta*ct*gamma^2*k + 12*Q*beta*ct*gamma^2 - 32*Q*beta*ct*gamma*k + 8*Q*beta*ct*gamma - 8*Q*beta*ct*k + 2*Q*beta*ct + 72*Q*beta*delta^6*w + 252*Q*beta*delta^5*gamma*w + 300*Q*beta*delta^5*w + 324*Q*beta*delta^4*gamma^2*w + 816*Q*beta*delta^4*gamma*w + 500*Q*beta*delta^4*w + 180*Q*beta*delta^3*gamma^3*w + 756*Q*beta*delta^3*gamma^2*w + 1000*Q*beta*delta^3*gamma*w + 424*Q*beta*delta^3*w + 36*Q*beta*delta^2*gamma^4*w + 264*Q*beta*delta^2*gamma^3*w + 612*Q*beta*delta^2*gamma^2*w + 576*Q*beta*delta^2*gamma*w + 192*Q*beta*delta^2*w + 24*Q*beta*delta*gamma^4*w + 116*Q*beta*delta*gamma^3*w + 204*Q*beta*delta*gamma^2*w + 156*Q*beta*delta*gamma*w + 44*Q*beta*delta*w + 4*Q*beta*gamma^4*w + 16*Q*beta*gamma^3*w + 24*Q*beta*gamma^2*w + 16*Q*beta*gamma*w + 4*Q*beta*w - 20*beta^2*c^2*delta^6 + 8*beta^2*c^2*delta^5*gamma - 80*beta^2*c^2*delta^5 + 96*beta^2*c^2*delta^4*gamma^2 + 40*beta^2*c^2*delta^4*gamma - 92*beta^2*c^2*delta^4 + 80*beta^2*c^2*delta^3*gamma^3 + 212*beta^2*c^2*delta^3*gamma^2 + 120*beta^2*c^2*delta^3*gamma - 12*beta^2*c^2*delta^3 + 16*beta^2*c^2*delta^2*gamma^4 + 112*beta^2*c^2*delta^2*gamma^3 + 216*beta^2*c^2*delta^2*gamma^2 + 160*beta^2*c^2*delta^2*gamma + 40*beta^2*c^2*delta^2 + 16*beta^2*c^2*delta*gamma^4 + 72*beta^2*c^2*delta*gamma^3 + 120*beta^2*c^2*delta*gamma^2 + 88*beta^2*c^2*delta*gamma + 24*beta^2*c^2*delta + 4*beta^2*c^2*gamma^4 + 16*beta^2*c^2*gamma^3 + 24*beta^2*c^2*gamma^2 + 16*beta^2*c^2*gamma + 4*beta^2*c^2 + 36*beta^2*c*co*delta^6 + 108*beta^2*c*co*delta^5*gamma + 164*beta^2*c*co*delta^5 + 108*beta^2*c*co*delta^4*gamma^2 + 408*beta^2*c*co*delta^4*gamma + 320*beta^2*c*co*delta^4 + 36*beta^2*c*co*delta^3*gamma^3 + 344*beta^2*c*co*delta^3*gamma^2 + 648*beta^2*c*co*delta^3*gamma + 340*beta^2*c*co*delta^3 + 120*beta^2*c*co*delta^2*gamma^3 + 444*beta^2*c*co*delta^2*gamma^2 + 528*beta^2*c*co*delta^2*gamma + 204*beta^2*c*co*delta^2 + 20*beta^2*c*co*delta*gamma^4 + 124*beta^2*c*co*delta*gamma^3 + 252*beta^2*c*co*delta*gamma^2 + 212*beta^2*c*co*delta*gamma + 64*beta^2*c*co*delta + 8*beta^2*c*co*gamma^4 + 32*beta^2*c*co*gamma^3 + 48*beta^2*c*co*gamma^2 + 32*beta^2*c*co*gamma + 8*beta^2*c*co + 36*beta^2*c*ct*delta^6 + 108*beta^2*c*ct*delta^5*gamma + 164*beta^2*c*ct*delta^5 + 108*beta^2*c*ct*delta^4*gamma^2 + 408*beta^2*c*ct*delta^4*gamma + 320*beta^2*c*ct*delta^4 + 36*beta^2*c*ct*delta^3*gamma^3 + 344*beta^2*c*ct*delta^3*gamma^2 + 648*beta^2*c*ct*delta^3*gamma + 340*beta^2*c*ct*delta^3 + 120*beta^2*c*ct*delta^2*gamma^3 + 444*beta^2*c*ct*delta^2*gamma^2 + 528*beta^2*c*ct*delta^2*gamma + 204*beta^2*c*ct*delta^2 + 20*beta^2*c*ct*delta*gamma^4 + 124*beta^2*c*ct*delta*gamma^3 + 252*beta^2*c*ct*delta*gamma^2 + 212*beta^2*c*ct*delta*gamma + 64*beta^2*c*ct*delta + 8*beta^2*c*ct*gamma^4 + 32*beta^2*c*ct*gamma^3 + 48*beta^2*c*ct*gamma^2 + 32*beta^2*c*ct*gamma + 8*beta^2*c*ct + 112*beta^2*c*delta^6*w + 200*beta^2*c*delta^5*gamma*w + 488*beta^2*c*delta^5*w + 24*beta^2*c*delta^4*gamma^2*w + 736*beta^2*c*delta^4*gamma*w + 824*beta^2*c*delta^4*w - 88*beta^2*c*delta^3*gamma^3*w + 264*beta^2*c*delta^3*gamma^2*w + 1056*beta^2*c*delta^3*gamma*w + 704*beta^2*c*delta^3*w - 32*beta^2*c*delta^2*gamma^4*w + 16*beta^2*c*delta^2*gamma^3*w + 456*beta^2*c*delta^2*gamma^2*w + 736*beta^2*c*delta^2*gamma*w + 328*beta^2*c*delta^2*w + 8*beta^2*c*delta*gamma^4*w + 104*beta^2*c*delta*gamma^3*w + 264*beta^2*c*delta*gamma^2*w + 248*beta^2*c*delta*gamma*w + 80*beta^2*c*delta*w + 8*beta^2*c*gamma^4*w + 32*beta^2*c*gamma^3*w + 48*beta^2*c*gamma^2*w + 32*beta^2*c*gamma*w + 8*beta^2*c*w - 9*beta^2*co^2*delta^6 - 36*beta^2*co^2*delta^5*gamma - 26*beta^2*co^2*delta^5 - 54*beta^2*co^2*delta^4*gamma^2 - 72*beta^2*co^2*delta^4*gamma - 11*beta^2*co^2*delta^4 - 36*beta^2*co^2*delta^3*gamma^3 - 62*beta^2*co^2*delta^3*gamma^2 + 8*beta^2*co^2*delta^3*gamma + 34*beta^2*co^2*delta^3 - 9*beta^2*co^2*delta^2*gamma^4 - 12*beta^2*co^2*delta^2*gamma^3 + 48*beta^2*co^2*delta^2*gamma^2 + 96*beta^2*co^2*delta^2*gamma + 45*beta^2*co^2*delta^2 + 4*beta^2*co^2*delta*gamma^4 + 32*beta^2*co^2*delta*gamma^3 + 72*beta^2*co^2*delta*gamma^2 + 64*beta^2*co^2*delta*gamma + 20*beta^2*co^2*delta + 3*beta^2*co^2*gamma^4 + 12*beta^2*co^2*gamma^3 + 18*beta^2*co^2*gamma^2 + 12*beta^2*co^2*gamma + 3*beta^2*co^2 - 18*beta^2*co*ct*delta^6 - 72*beta^2*co*ct*delta^5*gamma - 84*beta^2*co*ct*delta^5 - 108*beta^2*co*ct*delta^4*gamma^2 - 264*beta^2*co*ct*delta^4*gamma - 158*beta^2*co*ct*delta^4 - 72*beta^2*co*ct*delta^3*gamma^3 - 288*beta^2*co*ct*delta^3*gamma^2 - 368*beta^2*co*ct*delta^3*gamma - 152*beta^2*co*ct*delta^3 - 18*beta^2*co*ct*delta^2*gamma^4 - 120*beta^2*co*ct*delta^2*gamma^3 - 264*beta^2*co*ct*delta^2*gamma^2 - 240*beta^2*co*ct*delta^2*gamma - 78*beta^2*co*ct*delta^2 - 12*beta^2*co*ct*delta*gamma^4 - 56*beta^2*co*ct*delta*gamma^3 - 96*beta^2*co*ct*delta*gamma^2 - 72*beta^2*co*ct*delta*gamma - 20*beta^2*co*ct*delta - 2*beta^2*co*ct*gamma^4 - 8*beta^2*co*ct*gamma^3 - 12*beta^2*co*ct*gamma^2 - 8*beta^2*co*ct*gamma - 2*beta^2*co*ct - 72*beta^2*co*delta^6*w - 252*beta^2*co*delta^5*gamma*w - 300*beta^2*co*delta^5*w - 324*beta^2*co*delta^4*gamma^2*w - 816*beta^2*co*delta^4*gamma*w - 500*beta^2*co*delta^4*w - 180*beta^2*co*delta^3*gamma^3*w - 756*beta^2*co*delta^3*gamma^2*w - 1000*beta^2*co*delta^3*gamma*w - 424*beta^2*co*delta^3*w - 36*beta^2*co*delta^2*gamma^4*w - 264*beta^2*co*delta^2*gamma^3*w - 612*beta^2*co*delta^2*gamma^2*w - 576*beta^2*co*delta^2*gamma*w - 192*beta^2*co*delta^2*w - 24*beta^2*co*delta*gamma^4*w - 116*beta^2*co*delta*gamma^3*w - 204*beta^2*co*delta*gamma^2*w - 156*beta^2*co*delta*gamma*w - 44*beta^2*co*delta*w - 4*beta^2*co*gamma^4*w - 16*beta^2*co*gamma^3*w - 24*beta^2*co*gamma^2*w - 16*beta^2*co*gamma*w - 4*beta^2*co*w - 9*beta^2*ct^2*delta^6 - 36*beta^2*ct^2*delta^5*gamma - 26*beta^2*ct^2*delta^5 - 54*beta^2*ct^2*delta^4*gamma^2 - 72*beta^2*ct^2*delta^4*gamma - 11*beta^2*ct^2*delta^4 - 36*beta^2*ct^2*delta^3*gamma^3 - 62*beta^2*ct^2*delta^3*gamma^2 + 8*beta^2*ct^2*delta^3*gamma + 34*beta^2*ct^2*delta^3 - 9*beta^2*ct^2*delta^2*gamma^4 - 12*beta^2*ct^2*delta^2*gamma^3 + 48*beta^2*ct^2*delta^2*gamma^2 + 96*beta^2*ct^2*delta^2*gamma + 45*beta^2*ct^2*delta^2 + 4*beta^2*ct^2*delta*gamma^4 + 32*beta^2*ct^2*delta*gamma^3 + 72*beta^2*ct^2*delta*gamma^2 + 64*beta^2*ct^2*delta*gamma + 20*beta^2*ct^2*delta + 3*beta^2*ct^2*gamma^4 + 12*beta^2*ct^2*gamma^3 + 18*beta^2*ct^2*gamma^2 + 12*beta^2*ct^2*gamma + 3*beta^2*ct^2 - 72*beta^2*ct*delta^6*w - 252*beta^2*ct*delta^5*gamma*w - 300*beta^2*ct*delta^5*w - 324*beta^2*ct*delta^4*gamma^2*w - 816*beta^2*ct*delta^4*gamma*w - 500*beta^2*ct*delta^4*w - 180*beta^2*ct*delta^3*gamma^3*w - 756*beta^2*ct*delta^3*gamma^2*w - 1000*beta^2*ct*delta^3*gamma*w - 424*beta^2*ct*delta^3*w - 36*beta^2*ct*delta^2*gamma^4*w - 264*beta^2*ct*delta^2*gamma^3*w - 612*beta^2*ct*delta^2*gamma^2*w - 576*beta^2*ct*delta^2*gamma*w - 192*beta^2*ct*delta^2*w - 24*beta^2*ct*delta*gamma^4*w - 116*beta^2*ct*delta*gamma^3*w - 204*beta^2*ct*delta*gamma^2*w - 156*beta^2*ct*delta*gamma*w - 44*beta^2*ct*delta*w - 4*beta^2*ct*gamma^4*w - 16*beta^2*ct*gamma^3*w - 24*beta^2*ct*gamma^2*w - 16*beta^2*ct*gamma*w - 4*beta^2*ct*w - 128*beta^2*delta^6*w^2 - 352*beta^2*delta^5*gamma*w^2 - 544*beta^2*delta^5*w^2 - 336*beta^2*delta^4*gamma^2*w^2 - 1184*beta^2*delta^4*gamma*w^2 - 912*beta^2*delta^4*w^2 - 136*beta^2*delta^3*gamma^3*w^2 - 888*beta^2*delta^3*gamma^2*w^2 - 1528*beta^2*delta^3*gamma*w^2 - 776*beta^2*delta^3*w^2 - 20*beta^2*delta^2*gamma^4*w^2 - 272*beta^2*delta^2*gamma^3*w^2 - 840*beta^2*delta^2*gamma^2*w^2 - 944*beta^2*delta^2*gamma*w^2 - 356*beta^2*delta^2*w^2 - 28*beta^2*delta*gamma^4*w^2 - 168*beta^2*delta*gamma^3*w^2 - 336*beta^2*delta*gamma^2*w^2 - 280*beta^2*delta*gamma*w^2 - 84*beta^2*delta*w^2 - 8*beta^2*gamma^4*w^2 - 32*beta^2*gamma^3*w^2 - 48*beta^2*gamma^2*w^2 - 32*beta^2*gamma*w^2 - 8*beta^2*w^2)/(8*beta*(delta + gamma + 1)^2*(8*delta^3 + 14*delta^2*gamma + 18*delta^2 + 5*delta*gamma^2 + 16*delta*gamma + 11*delta + 2*gamma^2 + 4*gamma + 2));

fprintf('-----------------------------------------------------我是三号分割线-----------------------------------------------------\n');
wn=roundn(wn,-2);fprintf(2,['wn=',num2str(wn)]);
pon=roundn(pon,-2);fprintf(2,['\npon=',num2str(pon)]);
ptn=roundn(ptn,-2);fprintf(2,['\nptn=',num2str(ptn)]);
don=roundn(don,-2);fprintf(2,['\ndon=',num2str(don)]);
dtn=roundn(dtn,-2);fprintf(2,['\ndtn=',num2str(dtn)]);
Umn=roundn(Umn,-2);fprintf(2,['\nUmn=',num2str(Umn)]);
Uon=roundn(Uon,-2);fprintf(2,['\nUon=',num2str(Uon)]);
Utn=roundn(Utn,-2);fprintf(2,['\nUtn=',num2str(Utn)]);
Un=roundn(Un,-2);fprintf(2,['\nUn=',num2str(Un)]);
disp(' ');
%% 5.2 双渠道供应链的批发价格折扣与服务成本共担协调契约分析
%{
注意：如果代码运行缺少文本等如小图标题、图中数字0.09和0.49，并非是代码不全所致，解决办法是在MATLAB图象窗口手动插入。
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 图2（a）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear;
figure(1);
subplot(2,2,1);
set(0,'defaultfigurecolor','w');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.01
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.01;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));%把命令行窗口的结果直接粘贴至这里。
    PImx(i)=delta*((beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi + (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - delta*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi - 1) - (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - (ct*phi*(2*delta + 1)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta(i)*w + 4*beta*delta*theta(i)*w + 2*beta*gamma*theta(i)*w))/(4*delta + 2*gamma + 2);
end
plot(theta, PImx, 'r+', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.50
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.50;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));
    PImx(i)=delta*((beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi + (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - delta*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi - 1) - (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - (ct*phi*(2*delta + 1)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta(i)*w + 4*beta*delta*theta(i)*w + 2*beta*gamma*theta(i)*w))/(4*delta + 2*gamma + 2);
end
plot(theta, PImx, 'k.-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.99
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.99;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));
    PImx(i)=delta*((beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi + (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - delta*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi - 1) - (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - (ct*phi*(2*delta + 1)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta(i)*w + 4*beta*delta*theta(i)*w + 2*beta*gamma*theta(i)*w))/(4*delta + 2*gamma + 2);
end
plot(theta, PImx, 'b.', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PImn=729.19
plot([0,1],[729.19,729.19], 'g-', 'LineWidth', 1);
hold on;
axis([0 1 -2000 4000]);
set(gca,'XTick',[0:0.1:1],'YTick',[-2000:1000:4000]);
xlabel('折扣比例\theta','FontWeight','bold');ylabel('契约协调后制造商利润','FontWeight','bold');
legend({'服务成本共担比例为0.01','服务成本共担比例为0.50','服务成本共担比例为0.99','协同前制造商利润为729.19'});
legend boxoff;
text(0.504,-3417.721,-0.2000,'(a)制造商利润的变化趋势');
plot([0.15,0.15],[-2000,4000], 'g-', 'LineWidth', 1);
plot([0.21,0.21],[-2000,4000], 'g-', 'LineWidth', 1);
plot([0.27,0.27],[-2000,4000], 'g-', 'LineWidth', 1);
text(0.15,-1822.23,0.3,'$0.15$','interpreter','latex');
text(0.214,-1836.249,0.3,'$0.21$','interpreter','latex');
text(0.286,-1850.267,0.3,'$0.27$','interpreter','latex');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 图2（b）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear;
subplot(2,2,2);
set(0,'defaultfigurecolor','w');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.01
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.01;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));%把命令行窗口的结果直接粘贴至这里。
    PIox(i)=(ct*delta*phi*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta(i)*w + 4*beta*delta*theta(i)*w + 2*beta*gamma*theta(i)*w))/(4*delta + 2*gamma + 2) - gamma*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi - 1) - (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - ((beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(delta + gamma + 1)*(co + ct*phi + (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1)));
end
plot(theta, PIox, 'r+', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.50
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.50;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));
    PIox(i)=(ct*delta*phi*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta(i)*w + 4*beta*delta*theta(i)*w + 2*beta*gamma*theta(i)*w))/(4*delta + 2*gamma + 2) - gamma*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi - 1) - (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - ((beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(delta + gamma + 1)*(co + ct*phi + (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1)));
end
plot(theta, PIox, 'k.-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.99
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.99;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));
    PIox(i)=(ct*delta*phi*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta(i)*w + 4*beta*delta*theta(i)*w + 2*beta*gamma*theta(i)*w))/(4*delta + 2*gamma + 2) - gamma*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi - 1) - (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - ((beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(delta + gamma + 1)*(co + ct*phi + (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1)));
end
plot(theta, PIox, 'b.', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PIon=654.05
plot([0,1],[654.05,654.05], 'g-', 'LineWidth', 1);
hold on;
axis([0 1 600 2400]);
set(gca,'XTick',[0:0.1:1],'YTick',[600:200:2400]);
xlabel('折扣比例\theta','FontWeight','bold');ylabel('契约协调后线上渠道商利润','FontWeight','bold');
legend({'服务成本共担比例为0.01','服务成本共担比例为0.50','服务成本共担比例为0.99','协同前线上渠道商利润为654.05'});
legend boxoff;
text(1.678,-3341.772,-0.2000,'(b)线上渠道商利润的变化趋势');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 图2（c）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear;
subplot(2,2,3);
set(0,'defaultfigurecolor','w');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.01
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.01;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));%把命令行窗口的结果直接粘贴至这里。
    PItx(i)=(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(delta + gamma + 1)*(ct*(phi - 1) - (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + gamma*((beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi + (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + (ct*delta*phi*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta(i)*w + 4*beta*delta*theta(i)*w + 2*beta*gamma*theta(i)*w))/(4*delta + 2*gamma + 2);
end
plot(theta, PItx, 'r+', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.50
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.50;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));
    PItx(i)=(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(delta + gamma + 1)*(ct*(phi - 1) - (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + gamma*((beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi + (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + (ct*delta*phi*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta(i)*w + 4*beta*delta*theta(i)*w + 2*beta*gamma*theta(i)*w))/(4*delta + 2*gamma + 2);
end
plot(theta, PItx, 'k.-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.99
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.99;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));
    PItx(i)=(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(delta + gamma + 1)*(ct*(phi - 1) - (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi + theta(i)*w) - beta*delta*(c - theta(i)*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + gamma*((beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi + (c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(2*delta + gamma + 1) + (beta*delta*(c - theta(i)*w) - beta*(delta + gamma + 1)*(co + ct*phi + theta(i)*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + (ct*delta*phi*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta(i)*w + 4*beta*delta*theta(i)*w + 2*beta*gamma*theta(i)*w))/(4*delta + 2*gamma + 2);
end
plot(theta, PItx, 'b.', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PItn=-539.28
plot([0,1],[-539.28,-539.28], 'g-', 'LineWidth', 1);
hold on;
axis([0 1 -1000 500]);
set(gca,'XTick',[0:0.1:1],'YTick',[-1000:500:500]);
xlabel('折扣比例\theta','FontWeight','bold');ylabel('契约协调后线下渠道商利润','FontWeight','bold');
legend({'服务成本共担比例为0.01','服务成本共担比例为0.50','服务成本共担比例为0.99','协同前线下渠道商利润为-539.28'});
legend boxoff;
text(0.359 ,-11658.227 ,-0.2000,'(c)线下渠道商利润的变化趋势');
plot([0.09,0.09],[-1000,500], 'g-', 'LineWidth', 1);
text(0.15,-1822.23,0.3,'$0.09$','interpreter','latex');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 图2（d）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear;
subplot(2,2,4);
set(0,'defaultfigurecolor','w');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.01
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.01;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));%把命令行窗口的结果直接粘贴至这里。
    PIx(i)=(2*Q^2*delta^2*k^2 - 2*Q^2*delta^2*k + Q^2*delta^2 + 4*Q^2*delta*gamma*k^2 - 4*Q^2*delta*gamma*k + 2*Q^2*delta*gamma + 4*Q^2*delta*k^2 - 4*Q^2*delta*k + 2*Q^2*delta + 2*Q^2*gamma^2*k^2 - 2*Q^2*gamma^2*k + Q^2*gamma^2 + 4*Q^2*gamma*k^2 - 4*Q^2*gamma*k + 2*Q^2*gamma + 2*Q^2*k^2 - 2*Q^2*k + Q^2 - 2*Q*beta*c*delta^2 - 4*Q*beta*c*delta*gamma - 4*Q*beta*c*delta - 2*Q*beta*c*gamma^2 - 4*Q*beta*c*gamma - 2*Q*beta*c + 2*Q*beta*co*delta^2*k - 2*Q*beta*co*delta^2 + 4*Q*beta*co*delta*gamma*k - 4*Q*beta*co*delta*gamma + 4*Q*beta*co*delta*k - 4*Q*beta*co*delta + 2*Q*beta*co*gamma^2*k - 2*Q*beta*co*gamma^2 + 4*Q*beta*co*gamma*k - 4*Q*beta*co*gamma + 2*Q*beta*co*k - 2*Q*beta*co + 4*Q*beta*ct*delta^2*k*phi - 2*Q*beta*ct*delta^2*k - 2*Q*beta*ct*delta^2*phi + 8*Q*beta*ct*delta*gamma*k*phi - 4*Q*beta*ct*delta*gamma*k - 4*Q*beta*ct*delta*gamma*phi + 8*Q*beta*ct*delta*k*phi - 4*Q*beta*ct*delta*k - 4*Q*beta*ct*delta*phi + 4*Q*beta*ct*gamma^2*k*phi - 2*Q*beta*ct*gamma^2*k - 2*Q*beta*ct*gamma^2*phi + 8*Q*beta*ct*gamma*k*phi - 4*Q*beta*ct*gamma*k - 4*Q*beta*ct*gamma*phi + 4*Q*beta*ct*k*phi - 2*Q*beta*ct*k - 2*Q*beta*ct*phi - 6*beta^2*c^2*delta^2 - 4*beta^2*c^2*delta*gamma - 4*beta^2*c^2*delta + 2*beta^2*c*co*delta^2 + 4*beta^2*c*co*delta*gamma + 4*beta^2*c*co*delta + 2*beta^2*c*co*gamma^2 + 4*beta^2*c*co*gamma + 2*beta^2*c*co + 2*beta^2*c*ct*delta^2 + 4*beta^2*c*ct*delta*gamma + 4*beta^2*c*ct*delta + 2*beta^2*c*ct*gamma^2 + 4*beta^2*c*ct*gamma + 2*beta^2*c*ct + 16*beta^2*c*delta^2*theta(i)*w + 16*beta^2*c*delta*gamma*theta(i)*w + 16*beta^2*c*delta*theta(i)*w + 4*beta^2*c*gamma^2*theta(i)*w + 8*beta^2*c*gamma*theta(i)*w + 4*beta^2*c*theta(i)*w + beta^2*co^2*delta^2 + 2*beta^2*co^2*delta*gamma + 2*beta^2*co^2*delta + beta^2*co^2*gamma^2 + 2*beta^2*co^2*gamma + beta^2*co^2 + 2*beta^2*co*ct*delta^2*phi + 4*beta^2*co*ct*delta*gamma*phi + 4*beta^2*co*ct*delta*phi + 2*beta^2*co*ct*gamma^2*phi + 4*beta^2*co*ct*gamma*phi + 2*beta^2*co*ct*phi + 2*beta^2*ct^2*delta^2*phi^2 - 2*beta^2*ct^2*delta^2*phi + beta^2*ct^2*delta^2 + 4*beta^2*ct^2*delta*gamma*phi^2 - 4*beta^2*ct^2*delta*gamma*phi + 2*beta^2*ct^2*delta*gamma + 4*beta^2*ct^2*delta*phi^2 - 4*beta^2*ct^2*delta*phi + 2*beta^2*ct^2*delta + 2*beta^2*ct^2*gamma^2*phi^2 - 2*beta^2*ct^2*gamma^2*phi + beta^2*ct^2*gamma^2 + 4*beta^2*ct^2*gamma*phi^2 - 4*beta^2*ct^2*gamma*phi + 2*beta^2*ct^2*gamma + 2*beta^2*ct^2*phi^2 - 2*beta^2*ct^2*phi + beta^2*ct^2 - 8*beta^2*delta^2*theta(i)^2*w^2 - 8*beta^2*delta*gamma*theta(i)^2*w^2 - 8*beta^2*delta*theta(i)^2*w^2 - 2*beta^2*gamma^2*theta(i)^2*w^2 - 4*beta^2*gamma*theta(i)^2*w^2 - 2*beta^2*theta(i)^2*w^2)/(4*beta*(delta + gamma + 1)^2);
end
plot(theta, PIx, 'r+', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.50
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.50;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));
    PIx(i)=(2*Q^2*delta^2*k^2 - 2*Q^2*delta^2*k + Q^2*delta^2 + 4*Q^2*delta*gamma*k^2 - 4*Q^2*delta*gamma*k + 2*Q^2*delta*gamma + 4*Q^2*delta*k^2 - 4*Q^2*delta*k + 2*Q^2*delta + 2*Q^2*gamma^2*k^2 - 2*Q^2*gamma^2*k + Q^2*gamma^2 + 4*Q^2*gamma*k^2 - 4*Q^2*gamma*k + 2*Q^2*gamma + 2*Q^2*k^2 - 2*Q^2*k + Q^2 - 2*Q*beta*c*delta^2 - 4*Q*beta*c*delta*gamma - 4*Q*beta*c*delta - 2*Q*beta*c*gamma^2 - 4*Q*beta*c*gamma - 2*Q*beta*c + 2*Q*beta*co*delta^2*k - 2*Q*beta*co*delta^2 + 4*Q*beta*co*delta*gamma*k - 4*Q*beta*co*delta*gamma + 4*Q*beta*co*delta*k - 4*Q*beta*co*delta + 2*Q*beta*co*gamma^2*k - 2*Q*beta*co*gamma^2 + 4*Q*beta*co*gamma*k - 4*Q*beta*co*gamma + 2*Q*beta*co*k - 2*Q*beta*co + 4*Q*beta*ct*delta^2*k*phi - 2*Q*beta*ct*delta^2*k - 2*Q*beta*ct*delta^2*phi + 8*Q*beta*ct*delta*gamma*k*phi - 4*Q*beta*ct*delta*gamma*k - 4*Q*beta*ct*delta*gamma*phi + 8*Q*beta*ct*delta*k*phi - 4*Q*beta*ct*delta*k - 4*Q*beta*ct*delta*phi + 4*Q*beta*ct*gamma^2*k*phi - 2*Q*beta*ct*gamma^2*k - 2*Q*beta*ct*gamma^2*phi + 8*Q*beta*ct*gamma*k*phi - 4*Q*beta*ct*gamma*k - 4*Q*beta*ct*gamma*phi + 4*Q*beta*ct*k*phi - 2*Q*beta*ct*k - 2*Q*beta*ct*phi - 6*beta^2*c^2*delta^2 - 4*beta^2*c^2*delta*gamma - 4*beta^2*c^2*delta + 2*beta^2*c*co*delta^2 + 4*beta^2*c*co*delta*gamma + 4*beta^2*c*co*delta + 2*beta^2*c*co*gamma^2 + 4*beta^2*c*co*gamma + 2*beta^2*c*co + 2*beta^2*c*ct*delta^2 + 4*beta^2*c*ct*delta*gamma + 4*beta^2*c*ct*delta + 2*beta^2*c*ct*gamma^2 + 4*beta^2*c*ct*gamma + 2*beta^2*c*ct + 16*beta^2*c*delta^2*theta(i)*w + 16*beta^2*c*delta*gamma*theta(i)*w + 16*beta^2*c*delta*theta(i)*w + 4*beta^2*c*gamma^2*theta(i)*w + 8*beta^2*c*gamma*theta(i)*w + 4*beta^2*c*theta(i)*w + beta^2*co^2*delta^2 + 2*beta^2*co^2*delta*gamma + 2*beta^2*co^2*delta + beta^2*co^2*gamma^2 + 2*beta^2*co^2*gamma + beta^2*co^2 + 2*beta^2*co*ct*delta^2*phi + 4*beta^2*co*ct*delta*gamma*phi + 4*beta^2*co*ct*delta*phi + 2*beta^2*co*ct*gamma^2*phi + 4*beta^2*co*ct*gamma*phi + 2*beta^2*co*ct*phi + 2*beta^2*ct^2*delta^2*phi^2 - 2*beta^2*ct^2*delta^2*phi + beta^2*ct^2*delta^2 + 4*beta^2*ct^2*delta*gamma*phi^2 - 4*beta^2*ct^2*delta*gamma*phi + 2*beta^2*ct^2*delta*gamma + 4*beta^2*ct^2*delta*phi^2 - 4*beta^2*ct^2*delta*phi + 2*beta^2*ct^2*delta + 2*beta^2*ct^2*gamma^2*phi^2 - 2*beta^2*ct^2*gamma^2*phi + beta^2*ct^2*gamma^2 + 4*beta^2*ct^2*gamma*phi^2 - 4*beta^2*ct^2*gamma*phi + 2*beta^2*ct^2*gamma + 2*beta^2*ct^2*phi^2 - 2*beta^2*ct^2*phi + beta^2*ct^2 - 8*beta^2*delta^2*theta(i)^2*w^2 - 8*beta^2*delta*gamma*theta(i)^2*w^2 - 8*beta^2*delta*theta(i)^2*w^2 - 2*beta^2*gamma^2*theta(i)^2*w^2 - 4*beta^2*gamma*theta(i)^2*w^2 - 2*beta^2*theta(i)^2*w^2)/(4*beta*(delta + gamma + 1)^2);
end
plot(theta, PIx, 'k.-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi=0.99
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,phi=0.99;
theta = [0: 0.1: 1];
for i = 1: length(theta)
    w=(c + 2*c*delta + c*gamma + ct*phi + ct*delta*phi + ct*gamma*phi)/(theta(i)*(2*delta + gamma + 1));
    PIx(i)=(2*Q^2*delta^2*k^2 - 2*Q^2*delta^2*k + Q^2*delta^2 + 4*Q^2*delta*gamma*k^2 - 4*Q^2*delta*gamma*k + 2*Q^2*delta*gamma + 4*Q^2*delta*k^2 - 4*Q^2*delta*k + 2*Q^2*delta + 2*Q^2*gamma^2*k^2 - 2*Q^2*gamma^2*k + Q^2*gamma^2 + 4*Q^2*gamma*k^2 - 4*Q^2*gamma*k + 2*Q^2*gamma + 2*Q^2*k^2 - 2*Q^2*k + Q^2 - 2*Q*beta*c*delta^2 - 4*Q*beta*c*delta*gamma - 4*Q*beta*c*delta - 2*Q*beta*c*gamma^2 - 4*Q*beta*c*gamma - 2*Q*beta*c + 2*Q*beta*co*delta^2*k - 2*Q*beta*co*delta^2 + 4*Q*beta*co*delta*gamma*k - 4*Q*beta*co*delta*gamma + 4*Q*beta*co*delta*k - 4*Q*beta*co*delta + 2*Q*beta*co*gamma^2*k - 2*Q*beta*co*gamma^2 + 4*Q*beta*co*gamma*k - 4*Q*beta*co*gamma + 2*Q*beta*co*k - 2*Q*beta*co + 4*Q*beta*ct*delta^2*k*phi - 2*Q*beta*ct*delta^2*k - 2*Q*beta*ct*delta^2*phi + 8*Q*beta*ct*delta*gamma*k*phi - 4*Q*beta*ct*delta*gamma*k - 4*Q*beta*ct*delta*gamma*phi + 8*Q*beta*ct*delta*k*phi - 4*Q*beta*ct*delta*k - 4*Q*beta*ct*delta*phi + 4*Q*beta*ct*gamma^2*k*phi - 2*Q*beta*ct*gamma^2*k - 2*Q*beta*ct*gamma^2*phi + 8*Q*beta*ct*gamma*k*phi - 4*Q*beta*ct*gamma*k - 4*Q*beta*ct*gamma*phi + 4*Q*beta*ct*k*phi - 2*Q*beta*ct*k - 2*Q*beta*ct*phi - 6*beta^2*c^2*delta^2 - 4*beta^2*c^2*delta*gamma - 4*beta^2*c^2*delta + 2*beta^2*c*co*delta^2 + 4*beta^2*c*co*delta*gamma + 4*beta^2*c*co*delta + 2*beta^2*c*co*gamma^2 + 4*beta^2*c*co*gamma + 2*beta^2*c*co + 2*beta^2*c*ct*delta^2 + 4*beta^2*c*ct*delta*gamma + 4*beta^2*c*ct*delta + 2*beta^2*c*ct*gamma^2 + 4*beta^2*c*ct*gamma + 2*beta^2*c*ct + 16*beta^2*c*delta^2*theta(i)*w + 16*beta^2*c*delta*gamma*theta(i)*w + 16*beta^2*c*delta*theta(i)*w + 4*beta^2*c*gamma^2*theta(i)*w + 8*beta^2*c*gamma*theta(i)*w + 4*beta^2*c*theta(i)*w + beta^2*co^2*delta^2 + 2*beta^2*co^2*delta*gamma + 2*beta^2*co^2*delta + beta^2*co^2*gamma^2 + 2*beta^2*co^2*gamma + beta^2*co^2 + 2*beta^2*co*ct*delta^2*phi + 4*beta^2*co*ct*delta*gamma*phi + 4*beta^2*co*ct*delta*phi + 2*beta^2*co*ct*gamma^2*phi + 4*beta^2*co*ct*gamma*phi + 2*beta^2*co*ct*phi + 2*beta^2*ct^2*delta^2*phi^2 - 2*beta^2*ct^2*delta^2*phi + beta^2*ct^2*delta^2 + 4*beta^2*ct^2*delta*gamma*phi^2 - 4*beta^2*ct^2*delta*gamma*phi + 2*beta^2*ct^2*delta*gamma + 4*beta^2*ct^2*delta*phi^2 - 4*beta^2*ct^2*delta*phi + 2*beta^2*ct^2*delta + 2*beta^2*ct^2*gamma^2*phi^2 - 2*beta^2*ct^2*gamma^2*phi + beta^2*ct^2*gamma^2 + 4*beta^2*ct^2*gamma*phi^2 - 4*beta^2*ct^2*gamma*phi + 2*beta^2*ct^2*gamma + 2*beta^2*ct^2*phi^2 - 2*beta^2*ct^2*phi + beta^2*ct^2 - 8*beta^2*delta^2*theta(i)^2*w^2 - 8*beta^2*delta*gamma*theta(i)^2*w^2 - 8*beta^2*delta*theta(i)^2*w^2 - 2*beta^2*gamma^2*theta(i)^2*w^2 - 4*beta^2*gamma*theta(i)^2*w^2 - 2*beta^2*theta(i)^2*w^2)/(4*beta*(delta + gamma + 1)^2);
end
plot(theta, PIx, 'b.', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PIc=1236.14
plot([0,1],[1236.14,1236.14], 'g-', 'LineWidth', 1);
hold on;
axis([0 1 0 3500]);
set(gca,'XTick',[0:0.1:1],'YTick',[0:500:3500]);
xlabel('折扣比例\theta','FontWeight','bold');ylabel('契约协调后供应链总利润','FontWeight','bold');
legend({'服务成本共担比例为0.01','服务成本共担比例为0.50','服务成本共担比例为0.99','协同前供应链总利润为1236.14'});
legend boxoff;
text(0.527 ,-3398.734 ,-0.2,'(d)供应链总利润的变化趋势');
plot([0.49,0.49],[0,3500], 'g-', 'LineWidth', 1);
text(1.689 ,-11620.253 ,-0.2000,'$0.49$','interpreter','latex');
hold on;
title('图2 价格折扣契约对双渠道供应链各主体和总利润的影响','FontWeight','bold','position',[0 0 0]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 图3（a）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear;
figure(2);
subplot(2,2,1);
set(0,'defaultfigurecolor','w');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.09
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.09;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PImx(i)=delta*((beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi(i) + (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - delta*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi(i) - 1) - (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - (ct*phi(i)*(2*delta + 1)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta*w + 4*beta*delta*theta*w + 2*beta*gamma*theta*w))/(4*delta + 2*gamma + 2);
end
plot(phi, PImx, 'r+', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.12
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.12;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PImx(i)=delta*((beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi(i) + (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - delta*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi(i) - 1) - (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - (ct*phi(i)*(2*delta + 1)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta*w + 4*beta*delta*theta*w + 2*beta*gamma*theta*w))/(4*delta + 2*gamma + 2);
end
plot(phi, PImx, 'k.-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.15
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.15;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PImx(i)=delta*((beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi(i) + (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - delta*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi(i) - 1) - (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - (ct*phi(i)*(2*delta + 1)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta*w + 4*beta*delta*theta*w + 2*beta*gamma*theta*w))/(4*delta + 2*gamma + 2);
end
plot(phi, PImx, 'g-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PImn=729.19
plot([0,1],[729.19,729.19], 'g-', 'LineWidth', 1);
hold on;
axis([0 1 -2000 4000]);
set(gca,'XTick',[0:0.1:1],'YTick',[-2000:1000:4000]);
xlabel('服务成本共担比例\phi','FontWeight','bold');ylabel('契约协调后制造商利润','FontWeight','bold');
legend({'价格折扣比例为0.09','价格折扣比例为0.12','价格折扣比例为0.15','协同前制造商利润为729.19'});
legend boxoff;
text(0.504,-3417.721,-0.2000,'(a)制造商利润的变化趋势');
plot([0.03,0.03],[-2000,4000], 'g-', 'LineWidth', 1);
text(0.15,-1822.23,0.3,'$0.03$','interpreter','latex');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 图3（b）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear;
subplot(2,2,2);
set(0,'defaultfigurecolor','w');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.09
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.09;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PIox(i)=(ct*delta*phi(i)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta*w + 4*beta*delta*theta*w + 2*beta*gamma*theta*w))/(4*delta + 2*gamma + 2) - gamma*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi(i) - 1) - (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - ((beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(delta + gamma + 1)*(co + ct*phi(i) + (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1)));
end
plot(phi, PIox, 'r+', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.12
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.12;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PIox(i)=(ct*delta*phi(i)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta*w + 4*beta*delta*theta*w + 2*beta*gamma*theta*w))/(4*delta + 2*gamma + 2) - gamma*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi(i) - 1) - (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - ((beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(delta + gamma + 1)*(co + ct*phi(i) + (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1)));
end
plot(phi, PIox, 'k.-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.15
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.15;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PIox(i)=(ct*delta*phi(i)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta*w + 4*beta*delta*theta*w + 2*beta*gamma*theta*w))/(4*delta + 2*gamma + 2) - gamma*(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(ct*(phi(i) - 1) - (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) - ((beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(delta + gamma + 1)*(co + ct*phi(i) + (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1)));
end
plot(phi, PIox, 'g-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PIon=654.05
plot([0,1],[654.05,654.05], 'g-', 'LineWidth', 1);
hold on;
axis([0 1 600 2000]);
set(gca,'XTick',[0:0.1:1],'YTick',[600:200:2000]);
xlabel('服务成本共担比例\phi','FontWeight','bold');ylabel('契约协调后线上渠道商利润','FontWeight','bold');
legend({'价格折扣比例为0.09','价格折扣比例为0.12','价格折扣比例为0.15','协同前制造商利润为654.05'});
legend boxoff;
text(0.504,-3417.721,-0.2000,'(b)线上渠道商利润的变化趋势');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 图3（c）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear;
subplot(2,2,3);
set(0,'defaultfigurecolor','w');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.09
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.09;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PItx(i)=(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(delta + gamma + 1)*(ct*(phi(i) - 1) - (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + gamma*((beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi(i) + (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + (ct*delta*phi(i)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta*w + 4*beta*delta*theta*w + 2*beta*gamma*theta*w))/(4*delta + 2*gamma + 2);
end
plot(phi, PItx, 'r+', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.12
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.12;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PItx(i)=(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(delta + gamma + 1)*(ct*(phi(i) - 1) - (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + gamma*((beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi(i) + (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + (ct*delta*phi(i)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta*w + 4*beta*delta*theta*w + 2*beta*gamma*theta*w))/(4*delta + 2*gamma + 2);
end
plot(phi, PItx, 'k.-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.15
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.15;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PItx(i)=(Q*k - (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*(delta + gamma + 1)))*(delta + gamma + 1)*(ct*(phi(i) - 1) - (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*(delta + gamma + 1)*(ct - ct*phi(i) + theta*w) - beta*delta*(c - theta*w) + Q*k*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + gamma*((beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*(delta + gamma + 1)) - Q*(k - 1))*(co + ct*phi(i) + (c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(2*delta + gamma + 1) + (beta*delta*(c - theta*w) - beta*(delta + gamma + 1)*(co + ct*phi(i) + theta*w) + Q*(k - 1)*(delta + gamma + 1))/(2*beta*(delta + gamma + 1))) + (ct*delta*phi(i)*(beta*co - Q*delta - Q*gamma - Q + beta*ct - 2*beta*c*delta + beta*co*delta + beta*ct*delta + beta*co*gamma + beta*ct*gamma + 2*beta*theta*w + 4*beta*delta*theta*w + 2*beta*gamma*theta*w))/(4*delta + 2*gamma + 2);
end
plot(phi, PItx, 'g-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PItn=-539.28
plot([0,1],[-539.28,-539.28], 'g-', 'LineWidth', 1);
hold on;
axis([0 1 -600 200]);
set(gca,'XTick',[0:0.1:1],'YTick',[-600:100:200]);
xlabel('服务成本共担比例\phi','FontWeight','bold');ylabel('契约协调后线下渠道商利润','FontWeight','bold');
legend({'价格折扣比例为0.09','价格折扣比例为0.12','价格折扣比例为0.15','协同前线下渠道商利润为654.05'});
legend boxoff;
text(0.504,-3417.721,-0.2000,'(c)线下渠道商利润的变化趋势');
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 图3（d）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear;
subplot(2,2,4);
set(0,'defaultfigurecolor','w');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.09
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.09;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PIx(i)=(2*Q^2*delta^2*k^2 - 2*Q^2*delta^2*k + Q^2*delta^2 + 4*Q^2*delta*gamma*k^2 - 4*Q^2*delta*gamma*k + 2*Q^2*delta*gamma + 4*Q^2*delta*k^2 - 4*Q^2*delta*k + 2*Q^2*delta + 2*Q^2*gamma^2*k^2 - 2*Q^2*gamma^2*k + Q^2*gamma^2 + 4*Q^2*gamma*k^2 - 4*Q^2*gamma*k + 2*Q^2*gamma + 2*Q^2*k^2 - 2*Q^2*k + Q^2 - 2*Q*beta*c*delta^2 - 4*Q*beta*c*delta*gamma - 4*Q*beta*c*delta - 2*Q*beta*c*gamma^2 - 4*Q*beta*c*gamma - 2*Q*beta*c + 2*Q*beta*co*delta^2*k - 2*Q*beta*co*delta^2 + 4*Q*beta*co*delta*gamma*k - 4*Q*beta*co*delta*gamma + 4*Q*beta*co*delta*k - 4*Q*beta*co*delta + 2*Q*beta*co*gamma^2*k - 2*Q*beta*co*gamma^2 + 4*Q*beta*co*gamma*k - 4*Q*beta*co*gamma + 2*Q*beta*co*k - 2*Q*beta*co + 4*Q*beta*ct*delta^2*k*phi(i) - 2*Q*beta*ct*delta^2*k - 2*Q*beta*ct*delta^2*phi(i) + 8*Q*beta*ct*delta*gamma*k*phi(i) - 4*Q*beta*ct*delta*gamma*k - 4*Q*beta*ct*delta*gamma*phi(i) + 8*Q*beta*ct*delta*k*phi(i) - 4*Q*beta*ct*delta*k - 4*Q*beta*ct*delta*phi(i) + 4*Q*beta*ct*gamma^2*k*phi(i) - 2*Q*beta*ct*gamma^2*k - 2*Q*beta*ct*gamma^2*phi(i) + 8*Q*beta*ct*gamma*k*phi(i) - 4*Q*beta*ct*gamma*k - 4*Q*beta*ct*gamma*phi(i) + 4*Q*beta*ct*k*phi(i) - 2*Q*beta*ct*k - 2*Q*beta*ct*phi(i) - 6*beta^2*c^2*delta^2 - 4*beta^2*c^2*delta*gamma - 4*beta^2*c^2*delta + 2*beta^2*c*co*delta^2 + 4*beta^2*c*co*delta*gamma + 4*beta^2*c*co*delta + 2*beta^2*c*co*gamma^2 + 4*beta^2*c*co*gamma + 2*beta^2*c*co + 2*beta^2*c*ct*delta^2 + 4*beta^2*c*ct*delta*gamma + 4*beta^2*c*ct*delta + 2*beta^2*c*ct*gamma^2 + 4*beta^2*c*ct*gamma + 2*beta^2*c*ct + 16*beta^2*c*delta^2*theta*w + 16*beta^2*c*delta*gamma*theta*w + 16*beta^2*c*delta*theta*w + 4*beta^2*c*gamma^2*theta*w + 8*beta^2*c*gamma*theta*w + 4*beta^2*c*theta*w + beta^2*co^2*delta^2 + 2*beta^2*co^2*delta*gamma + 2*beta^2*co^2*delta + beta^2*co^2*gamma^2 + 2*beta^2*co^2*gamma + beta^2*co^2 + 2*beta^2*co*ct*delta^2*phi(i) + 4*beta^2*co*ct*delta*gamma*phi(i) + 4*beta^2*co*ct*delta*phi(i) + 2*beta^2*co*ct*gamma^2*phi(i) + 4*beta^2*co*ct*gamma*phi(i) + 2*beta^2*co*ct*phi(i) + 2*beta^2*ct^2*delta^2*phi(i)^2 - 2*beta^2*ct^2*delta^2*phi(i) + beta^2*ct^2*delta^2 + 4*beta^2*ct^2*delta*gamma*phi(i)^2 - 4*beta^2*ct^2*delta*gamma*phi(i) + 2*beta^2*ct^2*delta*gamma + 4*beta^2*ct^2*delta*phi(i)^2 - 4*beta^2*ct^2*delta*phi(i) + 2*beta^2*ct^2*delta + 2*beta^2*ct^2*gamma^2*phi(i)^2 - 2*beta^2*ct^2*gamma^2*phi(i) + beta^2*ct^2*gamma^2 + 4*beta^2*ct^2*gamma*phi(i)^2 - 4*beta^2*ct^2*gamma*phi(i) + 2*beta^2*ct^2*gamma + 2*beta^2*ct^2*phi(i)^2 - 2*beta^2*ct^2*phi(i) + beta^2*ct^2 - 8*beta^2*delta^2*theta^2*w^2 - 8*beta^2*delta*gamma*theta^2*w^2 - 8*beta^2*delta*theta^2*w^2 - 2*beta^2*gamma^2*theta^2*w^2 - 4*beta^2*gamma*theta^2*w^2 - 2*beta^2*theta^2*w^2)/(4*beta*(delta + gamma + 1)^2);
end
plot(phi, PIx, 'r+', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.12
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.12;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PIx(i)=(2*Q^2*delta^2*k^2 - 2*Q^2*delta^2*k + Q^2*delta^2 + 4*Q^2*delta*gamma*k^2 - 4*Q^2*delta*gamma*k + 2*Q^2*delta*gamma + 4*Q^2*delta*k^2 - 4*Q^2*delta*k + 2*Q^2*delta + 2*Q^2*gamma^2*k^2 - 2*Q^2*gamma^2*k + Q^2*gamma^2 + 4*Q^2*gamma*k^2 - 4*Q^2*gamma*k + 2*Q^2*gamma + 2*Q^2*k^2 - 2*Q^2*k + Q^2 - 2*Q*beta*c*delta^2 - 4*Q*beta*c*delta*gamma - 4*Q*beta*c*delta - 2*Q*beta*c*gamma^2 - 4*Q*beta*c*gamma - 2*Q*beta*c + 2*Q*beta*co*delta^2*k - 2*Q*beta*co*delta^2 + 4*Q*beta*co*delta*gamma*k - 4*Q*beta*co*delta*gamma + 4*Q*beta*co*delta*k - 4*Q*beta*co*delta + 2*Q*beta*co*gamma^2*k - 2*Q*beta*co*gamma^2 + 4*Q*beta*co*gamma*k - 4*Q*beta*co*gamma + 2*Q*beta*co*k - 2*Q*beta*co + 4*Q*beta*ct*delta^2*k*phi(i) - 2*Q*beta*ct*delta^2*k - 2*Q*beta*ct*delta^2*phi(i) + 8*Q*beta*ct*delta*gamma*k*phi(i) - 4*Q*beta*ct*delta*gamma*k - 4*Q*beta*ct*delta*gamma*phi(i) + 8*Q*beta*ct*delta*k*phi(i) - 4*Q*beta*ct*delta*k - 4*Q*beta*ct*delta*phi(i) + 4*Q*beta*ct*gamma^2*k*phi(i) - 2*Q*beta*ct*gamma^2*k - 2*Q*beta*ct*gamma^2*phi(i) + 8*Q*beta*ct*gamma*k*phi(i) - 4*Q*beta*ct*gamma*k - 4*Q*beta*ct*gamma*phi(i) + 4*Q*beta*ct*k*phi(i) - 2*Q*beta*ct*k - 2*Q*beta*ct*phi(i) - 6*beta^2*c^2*delta^2 - 4*beta^2*c^2*delta*gamma - 4*beta^2*c^2*delta + 2*beta^2*c*co*delta^2 + 4*beta^2*c*co*delta*gamma + 4*beta^2*c*co*delta + 2*beta^2*c*co*gamma^2 + 4*beta^2*c*co*gamma + 2*beta^2*c*co + 2*beta^2*c*ct*delta^2 + 4*beta^2*c*ct*delta*gamma + 4*beta^2*c*ct*delta + 2*beta^2*c*ct*gamma^2 + 4*beta^2*c*ct*gamma + 2*beta^2*c*ct + 16*beta^2*c*delta^2*theta*w + 16*beta^2*c*delta*gamma*theta*w + 16*beta^2*c*delta*theta*w + 4*beta^2*c*gamma^2*theta*w + 8*beta^2*c*gamma*theta*w + 4*beta^2*c*theta*w + beta^2*co^2*delta^2 + 2*beta^2*co^2*delta*gamma + 2*beta^2*co^2*delta + beta^2*co^2*gamma^2 + 2*beta^2*co^2*gamma + beta^2*co^2 + 2*beta^2*co*ct*delta^2*phi(i) + 4*beta^2*co*ct*delta*gamma*phi(i) + 4*beta^2*co*ct*delta*phi(i) + 2*beta^2*co*ct*gamma^2*phi(i) + 4*beta^2*co*ct*gamma*phi(i) + 2*beta^2*co*ct*phi(i) + 2*beta^2*ct^2*delta^2*phi(i)^2 - 2*beta^2*ct^2*delta^2*phi(i) + beta^2*ct^2*delta^2 + 4*beta^2*ct^2*delta*gamma*phi(i)^2 - 4*beta^2*ct^2*delta*gamma*phi(i) + 2*beta^2*ct^2*delta*gamma + 4*beta^2*ct^2*delta*phi(i)^2 - 4*beta^2*ct^2*delta*phi(i) + 2*beta^2*ct^2*delta + 2*beta^2*ct^2*gamma^2*phi(i)^2 - 2*beta^2*ct^2*gamma^2*phi(i) + beta^2*ct^2*gamma^2 + 4*beta^2*ct^2*gamma*phi(i)^2 - 4*beta^2*ct^2*gamma*phi(i) + 2*beta^2*ct^2*gamma + 2*beta^2*ct^2*phi(i)^2 - 2*beta^2*ct^2*phi(i) + beta^2*ct^2 - 8*beta^2*delta^2*theta^2*w^2 - 8*beta^2*delta*gamma*theta^2*w^2 - 8*beta^2*delta*theta^2*w^2 - 2*beta^2*gamma^2*theta^2*w^2 - 4*beta^2*gamma*theta^2*w^2 - 2*beta^2*theta^2*w^2)/(4*beta*(delta + gamma + 1)^2);
end
plot(phi, PIx, 'k.-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=0.15
beta=0.7,c=10,Q=100,k=0.4,co=3,ct=5,delta=0.8,gamma=0.8,theta=0.15;
phi = [0: 0.1: 1];
for i = 1: length(phi)
    w=(c + 2*c*delta + c*gamma + ct*phi(i) + ct*delta*phi(i) + ct*gamma*phi(i))/(theta*(2*delta + gamma + 1));
    PIx(i)=(2*Q^2*delta^2*k^2 - 2*Q^2*delta^2*k + Q^2*delta^2 + 4*Q^2*delta*gamma*k^2 - 4*Q^2*delta*gamma*k + 2*Q^2*delta*gamma + 4*Q^2*delta*k^2 - 4*Q^2*delta*k + 2*Q^2*delta + 2*Q^2*gamma^2*k^2 - 2*Q^2*gamma^2*k + Q^2*gamma^2 + 4*Q^2*gamma*k^2 - 4*Q^2*gamma*k + 2*Q^2*gamma + 2*Q^2*k^2 - 2*Q^2*k + Q^2 - 2*Q*beta*c*delta^2 - 4*Q*beta*c*delta*gamma - 4*Q*beta*c*delta - 2*Q*beta*c*gamma^2 - 4*Q*beta*c*gamma - 2*Q*beta*c + 2*Q*beta*co*delta^2*k - 2*Q*beta*co*delta^2 + 4*Q*beta*co*delta*gamma*k - 4*Q*beta*co*delta*gamma + 4*Q*beta*co*delta*k - 4*Q*beta*co*delta + 2*Q*beta*co*gamma^2*k - 2*Q*beta*co*gamma^2 + 4*Q*beta*co*gamma*k - 4*Q*beta*co*gamma + 2*Q*beta*co*k - 2*Q*beta*co + 4*Q*beta*ct*delta^2*k*phi(i) - 2*Q*beta*ct*delta^2*k - 2*Q*beta*ct*delta^2*phi(i) + 8*Q*beta*ct*delta*gamma*k*phi(i) - 4*Q*beta*ct*delta*gamma*k - 4*Q*beta*ct*delta*gamma*phi(i) + 8*Q*beta*ct*delta*k*phi(i) - 4*Q*beta*ct*delta*k - 4*Q*beta*ct*delta*phi(i) + 4*Q*beta*ct*gamma^2*k*phi(i) - 2*Q*beta*ct*gamma^2*k - 2*Q*beta*ct*gamma^2*phi(i) + 8*Q*beta*ct*gamma*k*phi(i) - 4*Q*beta*ct*gamma*k - 4*Q*beta*ct*gamma*phi(i) + 4*Q*beta*ct*k*phi(i) - 2*Q*beta*ct*k - 2*Q*beta*ct*phi(i) - 6*beta^2*c^2*delta^2 - 4*beta^2*c^2*delta*gamma - 4*beta^2*c^2*delta + 2*beta^2*c*co*delta^2 + 4*beta^2*c*co*delta*gamma + 4*beta^2*c*co*delta + 2*beta^2*c*co*gamma^2 + 4*beta^2*c*co*gamma + 2*beta^2*c*co + 2*beta^2*c*ct*delta^2 + 4*beta^2*c*ct*delta*gamma + 4*beta^2*c*ct*delta + 2*beta^2*c*ct*gamma^2 + 4*beta^2*c*ct*gamma + 2*beta^2*c*ct + 16*beta^2*c*delta^2*theta*w + 16*beta^2*c*delta*gamma*theta*w + 16*beta^2*c*delta*theta*w + 4*beta^2*c*gamma^2*theta*w + 8*beta^2*c*gamma*theta*w + 4*beta^2*c*theta*w + beta^2*co^2*delta^2 + 2*beta^2*co^2*delta*gamma + 2*beta^2*co^2*delta + beta^2*co^2*gamma^2 + 2*beta^2*co^2*gamma + beta^2*co^2 + 2*beta^2*co*ct*delta^2*phi(i) + 4*beta^2*co*ct*delta*gamma*phi(i) + 4*beta^2*co*ct*delta*phi(i) + 2*beta^2*co*ct*gamma^2*phi(i) + 4*beta^2*co*ct*gamma*phi(i) + 2*beta^2*co*ct*phi(i) + 2*beta^2*ct^2*delta^2*phi(i)^2 - 2*beta^2*ct^2*delta^2*phi(i) + beta^2*ct^2*delta^2 + 4*beta^2*ct^2*delta*gamma*phi(i)^2 - 4*beta^2*ct^2*delta*gamma*phi(i) + 2*beta^2*ct^2*delta*gamma + 4*beta^2*ct^2*delta*phi(i)^2 - 4*beta^2*ct^2*delta*phi(i) + 2*beta^2*ct^2*delta + 2*beta^2*ct^2*gamma^2*phi(i)^2 - 2*beta^2*ct^2*gamma^2*phi(i) + beta^2*ct^2*gamma^2 + 4*beta^2*ct^2*gamma*phi(i)^2 - 4*beta^2*ct^2*gamma*phi(i) + 2*beta^2*ct^2*gamma + 2*beta^2*ct^2*phi(i)^2 - 2*beta^2*ct^2*phi(i) + beta^2*ct^2 - 8*beta^2*delta^2*theta^2*w^2 - 8*beta^2*delta*gamma*theta^2*w^2 - 8*beta^2*delta*theta^2*w^2 - 2*beta^2*gamma^2*theta^2*w^2 - 4*beta^2*gamma*theta^2*w^2 - 2*beta^2*theta^2*w^2)/(4*beta*(delta + gamma + 1)^2);
end
plot(phi, PIx, 'g-', 'LineWidth', 1);
hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PIc=1236.14
plot([0,1],[1236.14,1236.14], 'g-', 'LineWidth', 1);
hold on;
axis([0 1 1200 2800]);
set(gca,'XTick',[0:0.1:1],'YTick',[1200:200:2800]);
xlabel('服务成本共担比例\phi','FontWeight','bold');ylabel('契约协调后供应链总利润','FontWeight','bold');
legend({'价格折扣比例为0.09','价格折扣比例为0.12','价格折扣比例为0.15','协同前供应链总利润为1236.14'});
legend boxoff;
text(0.504,-3417.721,-0.2000,'(d)供应链总利润的变化趋势');
hold on;
title('图3 服务成本共担契约对双渠道供应链各主体和总利润的影响','FontWeight','bold','position',[0 0 0]);