ClearAll;
Needs["Statistics`DataManipulation`"];
Needs["Graphics`Graphics`"];
pairs = {{0,0},{1,167/256},{2,215/256},{3,238/256},{4,244/256},
		  {5,247/256},{6,253/256},{7,254/256},{8,255/256},{9,256/256}};
times = Flatten[pairs][[Table[2 i - 1,{i,10}]]];
deaths = Flatten[pairs][[Table[2 i ,{i,10}]]];

GGL[x_,a_,b_]:= a x + b x^2;

fl[a_,b_]:= Apply[Plus,(Table[GGL[i,a,b],{i,1,2}] - Drop[Drop[deaths,1],-7])^2];
fml:= FindMinimum[fl[x,y],{x,1},{y,1}];
a0 = x/.fml[[2]]; b0 = y/.fml[[2]];
GL[x_]:=  Evaluate[a0 x + b0 x^2];

GGH[x_,c_]:= 1 - c x^(-2);

fh[c_]:= Apply[Plus,(Table[GGH[i,c],{i,2,9}] - Drop[Drop[deaths,1],1])^2];
fmh:= FindMinimum[fh[x],{x,1}];
c0 = x/.fmh[[2]];

GH[x_]:=  Evaluate[1 - c0 x^(-2)];

xcross = xx/.FindRoot[GH[xx]==GL[xx],{xx,2,2.1}];

g01star[x_]:=  a0 + 2 b0 x; 
g02star[x_]:=  2 c0 x^(-3);

P1a:= Plot[GL[x],{x,0,xcross}];
P1b:= Plot[GH[x],{x,xcross,9}];

P1:= Show[P1a,P1b,PlotRange->{{0,9},{0,1}}];
P3:= Show[Plot[g01star[x],{x,0,xcross}],
		Plot[g02star[x],{x,xcross,9}],
		AxesOrigin -> {0,-0.1},
		PlotRange->{{0,9},{-0.1,1}},PlotLabel -> "(a)"];

G1[x_,d_,dd_]:= d x + ((1 -3 d - 20 dd)x^2)/8 + dd x^3;
G2[x_,d_,dd_]:= 1 - 2(1-d+4dd) x^(-2);

GG[x_,d_,dd_]:= Which[x > 2, G2[x,d,dd], x <= 2, G1[x,d,dd]];

f[d_,dd_]:= Apply[Plus,(Table[GG[i,d,dd],{i,1,9}] - Drop[deaths,1])^2];
fm:= FindMinimum[f[x,y],{x,1},{y,1}];
d0 = x/.fm[[2]]; dd0 = y/.fm[[2]];
g1star[x_]:=  d0 + (1 -3 d0 -20 dd0)x/4 + 3 dd0 x^2; 
g2star[x_]:=  4(1 - d0 + 4 dd0) x^(-3);
gstar[x_]:= Which[x > 2, g2star[x], x <= 2, g1star[x]];

G1star[x_]:=  Evaluate[G1[x,d0,dd0]]; 
G2star[x_]:=  Evaluate[G2[x,d0,dd0]];
Gstar[x_]:= Which[x > 2, G2star[x], x <= 2, G1star[x]];

Print[G1star[x],"    ",G2star[x],"    ",fm[[1]]];
upplot = Plot[Gstar[x],{x,0,9},PlotPoints->70,
	PlotRange->{{0,9.1},{0,1.05}},
	AxesLabel ->{"t","y"},
	PlotLabel->FontForm["y = W(t)",{"Courier",10}]];
up = Show[upplot,ListPlot[pairs]];
downplot = Plot[gstar[x],{x,0,9},PlotPoints->70,
	PlotRange->{{0,9.1},{0,1.2}},
	AxesLabel ->{"t","y"},
	PlotLabel->FontForm["y = W'(t)",{"Courier",10}]];
down = Show[downplot];
Print["For t <= 2 we have W'(t) = ",g1star[t]];
Print["For t >= 2 we have W'(t) = ",g2star[t]];
Print["So"];Print[" "];
Print["W(2-) = ",G1star[2]];Print["W(2+) = ",G2star[2]];
Print["W'(2-) = ",g1star[2]];Print["W'(2+) = ",g2star[2]];
Print[" "];Print["implying that W is smooth."];Print[" "];
Print[" "];
Print["For t <= 2 we have W''(t) = ",D[g1star[t],t]];
Print["For t >= 2 we have W''(t) = ",D[g2star[t],t]];
Print["So"];Print[" "];
Print["W''(2-) = ",D[g1star[t],t]/.t->2];
Print["W''(2+) = ",D[g2star[t],t]/.t->2];
Print[" "];Print["implying an inflection point where t = 2."];Print[" "];
Print["The other inflection point is where t = ",t/.Part[Solve[D[g1star[t],t]==0,t],1]];
Print[" "];Print["The answer to (iii) is No! Just look at these graphs!"];
Show[GraphicsArray[{{up},{down}}]]