% Read function definitions
prob9;

% Table 9.1
mth=2;   % Integrate using quad
[alf,f]=maximize(mth,7);
for k=1:7, disp(sprintf('%24.17f   %24.17f',alf(k),f(k))); end
disp('   Table 9.1\n')

% Table 9.2
[alf,f]=findzero(mth);
for k=1:length(alf), disp(sprintf('%24.16f   %24.16f',alf(k),f(k))); end
disp('   Table 9.2\n')

% Table 9.3
alf(2)=pi/4; alf(1)=alf(2)-0.005; global method bigM; 
method=3; bigM=16;
f(1)=dI(alf(1)); f(2)=dI(alf(2)); 
for k=2:10, 
  alf(k+1)=alf(k)-(alf(k)-alf(k-1))*f(k)/(f(k)-f(k-1));
  f(k+1)=dI(alf(k+1));
  if abs(f(k+1))>=abs(f(k)), break, end
end
for j=1:k+1, disp(sprintf('%24.16f   %24.16f',alf(j),f(j))); end
disp('   Table 9.3\n')