ClearAll;
$DefaultFont = {"Times-Roman",10};
Needs["Graphics`Graphics`"];
cumulative = {312/545,413/545,487/545,522/545,537/545,1};
datapoints = {{0,0},{1,312/545},{2,413/545},
			  {3,487/545},{4,522/545},{5,537/545},{6,1}};
density = {312/545,101/545,74/545,35/545,15/545,8/545};
F1[A_,th_,cc_,t_] = (3 th + (th + 2)A)t/(2(th+1)cc) - A(t/cc)^2 + ((-th + (th+2)A)/(2(th+3))) (t/cc)^3;
F2[A_,th_,cc_,t_] = 1 - ((3-A)/((3+th)(1+th))) (cc/t)^th;

FF[A_,th_,cc_,t_] = Which[t<=2,F1[A,th,cc,t],t>2,F2[A,th,cc,t]];

error[A_,th_,cc_] = Apply[Plus,Table[(FF[A,th,cc,n]-cumulative[[n]])^2,{n,1,6}]];
leasterror = FindMinimum[error[A,th,cc],{A,2,2.01,0,4},{th,2,2.01,0,4},{cc,2,2.01,0,4}];
c = cc/.leasterror[[2]]; K1 = 10; K2 = 20; Kmax = 30; a1 = 0.5; a2 = 2;
poly[a_,b_,f_,njags_]:= Join[Table[{a+i(b-a)/njags,f[a+i(b-a)/njags]},
				  {i,0,njags}],{{b,0},{a,0}}];
shade[a_,b_,c_,f_,njags_]:= Graphics[{GrayLevel[c],Polygon[poly[a,b,f,njags]]}];
fun1[z_]:= z^(-a1);
fun2[z_]:= z^(-a2);
panel1[K_]:= Show[shade[c,K,0.65,fun1,200],Plot[fun1[z],{z,c,K},PlotPoints->100],
				Axes->True, AxesOrigin->{-0.01Kmax,-0.025fun1[c]},
				Ticks->{{{c,"c"},{K,"K"}},{0,{fun1[c],"u(c)"}}},
				PlotRange->{{0,Kmax},{0,fun1[c]}},
				PlotLabel->StringForm["a = ``, K = ``",a1,K],
				AxesLabel->{"   t","y"}];
panel2[K_]:= Show[shade[c,K,0.65,fun2,200],Plot[fun2[z],{z,c,K},PlotPoints->100],
				Axes->True, AxesOrigin->{-0.01Kmax,-0.025fun2[c]},
				Ticks->{{{c,"c"},{K,"K"}},{0,{fun2[c],"u(c)"}}},
				PlotRange->{{0,Kmax},{0,fun2[c]}},
				PlotLabel->StringForm["a = ``, K = ``",a2,K],
				AxesLabel->{"   t","y"}];
Show[GraphicsArray[{{panel2[K1],panel1[K1]},{panel2[K2],panel1[K2]},{panel2[Kmax],panel1[Kmax]}}],
	 PlotLabel->"           Figure 27.3\n               y = u(t)\nShaded area = Int(u,[c,K])"]