Clear["Global`*"]
ClearAll[Subscript]

(*Model primitives*)
F[z_]:=z;
f[z_]:=D[F[z],z];
(UnderBar[ϵ])=0;
(OverBar[ϵ])=1;

(*Distribution of the second good quality Y conditional on observed quality X of the first one in the region Y>X*)
g[y_,x_]:=(((1-p)(1-F[x]))/((1-p)(1-F[x])+p F[x]))((f[y])/((1-F[x])))/.p->((1+ξ)/2);

(*Probability of visiting the second firm*)
Pr[z_]:=((1+ξ)/2)+((1-ξ)/2)(1-(((1-F[z]))^2));

(*Incremental utility of one extra search*)
ΔU[z_]:=FullSimplify[(∫_z^(OverBar[ϵ]))(y-z)g[y,z]\[DifferentialD]y];

(*Product Quality and Search Expenditures*)
MQ[z_] :=Pr[z](∫_(UnderBar[ϵ])^(OverBar[ϵ]))ϵ 2f[ϵ]F[ϵ]\[DifferentialD]ϵ+ (1-Pr[z])(∫_(UnderBar[ϵ])^(OverBar[ϵ]))ϵ 2f[ϵ](1-F[ϵ])\[DifferentialD]ϵ;
SE[z_]:=s(1+((1+ξ)/2)((F[z])^2)+((1-ξ)/2)(1-(((1-F[z]))^2)));

(*Demand and Price*)
Demand[p_,P_,z_]:=(1/2)(1+F[z])(1-F[z+p-P])+(∫_(UnderBar[ϵ])^(z+p-P))f[ϵ]F[ϵ-p+P]\[DifferentialD]ϵ;
Price[z_]:=p/.Solve[p==FullSimplify[-((Demand[p,p,z])/(D[Demand[p,P,z],p]/.P->p))],p][[1]];

(*Derivatives of main expressions*)
DPr[z_]:=FullSimplify[D[Pr[z[ξ]],ξ]/.z[ξ]->z];
DMQ[z_]:=DPr[z]((∫_(UnderBar[ϵ])^(OverBar[ϵ]))ϵ 2f[ϵ]F[ϵ]\[DifferentialD]ϵ-(∫_(UnderBar[ϵ])^(OverBar[ϵ]))ϵ 2f[ϵ](1-F[ϵ])\[DifferentialD]ϵ);
DSE[z_]:=Simplify[D[SE[z[ξ]],ξ]/.z[ξ]->z];
DTS[z_]:=DMQ[z]-DSE[z];

DPrice[z_]:=Simplify[D[Price[z[ξ]],ξ]/.z[ξ]->z];
DCS[z_]:=DTS[z]-DPrice[z];

(*Computation of X and ((∂X)/(∂ξ))*)
X=Normal[z/.Solve[s==ΔU[z],z,Reals][[1]]];
ϕ'[ξ]=FullSimplify[z'[ξ]/.Solve[FullSimplify[0==D[ΔU[z[ξ]],ξ]],z'[ξ]][[1]]/.z[ξ]->ϕ];

(*Graphs*)
(*Contourplots DMQ=0, DCS=0*)
FullSimplify[Simplify[DMQ[ϕ]/.ϕ->X]];
plot1=ContourPlot[{%==0},{ξ,0,1},{s,0,(1/2)},PlotPoints->200,FrameLabel->{"ξ",Rotate["s",-90Degree]},BoundaryStyle->{Black,Thick},PlotRange->{{0,1},{0,0.6}},ContourStyle->{{Black,Dashed}},FrameTicks->{{0,0.25,0.5,0.75,1},{.1,.2,.3,.4,.5,.6}},FrameTicksStyle->{{Black,White},{Black,White}}];
FullSimplify[Simplify[DCS[ϕ]/.ϕ->X]];
plot2=ContourPlot[{%==0},{ξ,0,1},{s,0,(1/2)},PlotPoints->100,FrameLabel->{"ξ",Rotate["s",-90Degree]},BoundaryStyle->{Black,Thick},PlotRange->{{0,1},{0,0.6}},ContourStyle->{{Black,Dotted}}];
plot3=Plot[{(1/2)},{ξ,0,1},PlotRange->{{0,1},{0,0.6}},PlotStyle->{{Gray,Thin}},Filling->Top];
text=Graphics[{Text["Product Quality ↓ \n Consumer Welfare↓",{0.15,0.035}],
Text["Product Quality ↑ \n Consumer Welfare↓",{0.6,0.17}],
Text["Product Quality ↑ \n Consumer Welfare↑",{0.4,0.4}],
Text["No Search \nProduct Quality ↑ \n Consumer Surplus↑",{0.5,0.55}]}];
Show[plot1,plot2,plot3(*,plot4*),text]

(*Contourplot DTS=0*)
s/.Normal[Solve[(DTS[ϕ]/.ϕ->X)==0&&0<ξ<1&&0<s<(1/2),s]][[1]];
plot21=Plot[%,{ξ,0,1},PlotRange->{{0,1},{0,0.02}},PlotStyle->{{Dashed,Black}},Axes->False,Frame->True,FrameLabel->{"ξ",Rotate["s",-90Degree]},PlotRange->{{0,1},{0,0.6}},FrameTicks->{{0,0.25,0.5,0.75,1},{0.01,0.02}},FrameTicksStyle->{{Black,White},{Black,White}},AspectRatio->1];
text2=Graphics[{Text["Total Welfare ↓",{0.15,0.0025}],Text["Total Welfare ↑",{0.5,0.01}]}];
Show[plot21,text2]

g[y_,x_]:=(((1-p)(1-F[x]))/((1-p)(1-F[x])+p F[x]))((f[y])/((1-F[x])))/.p->((1+q)/2);
test=((FullSimplify[g[ϵ,x]])/(f[ϵ]))
Plot3D[test,{q,0,1},{F[x],0,1},AxesLabel->Automatic]
ContourPlot[test,{q,0,1},{F[x],0,1},FrameLabel->Automatic,Contours->{0.99,0.95,.9,.8,.7,.6,.5,.4,.3,.2,.1},ContourLabels->True ]