function hw01p2()
% P.G. Martinsson, Jan 14, 2011.

a          = 0;
b          = 1;
n_quad_vec = 50.*(2.^(1:6));

f = inline('exp(t.*t).*cos(37*t)','t');
%f = inline('exp(t)','t');
%f = inline('cos(1./(0.01 + t.*t))','t');

t1_vec = zeros(size(n_quad_vec));
t2_vec = zeros(size(n_quad_vec));

for icount = 1:length(n_quad_vec)
  n_quad = n_quad_vec(icount);
  tic
  I1_vec(icount) = trapezoidal1(f,a,b,n_quad);
  t1_vec(icount) = toc;
  tic
  for j = 1:1000
    I2_vec(icount) = trapezoidal2(f,a,b,n_quad);
  end
  t2_vec(icount) = toc/1000;
end

figure(1)
loglog(n_quad_vec,t1_vec,'x-',...
       n_quad_vec,t2_vec,'x-');
title('Execution time for two versions of the trapezoidal rule')
legend('trapezoidal1','trapezoidal2')
xlabel('n')
ylabel('t')
axis tight
%print -depsc hw01p2_fig1.eps

figure(2)
loglog(n_quad_vec,t1_vec./n_quad_vec,'x-',...
       n_quad_vec,t2_vec./n_quad_vec,'x-')
title('Execution time for two versions of the trapezoidal rule')
legend('trapezoidal1','trapezoidal2')
xlabel('n')
ylabel('t/n')
axis tight
%print -depsc hw01p2_fig2.eps

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function I = trapezoidal1(f,a,b,n)
h = (b-a)/n;
I = 0.5*h*(f(a) + f(b));
for icount = 1:(n-1)
   I = I + h*f(a+icount*h);
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function I = trapezoidal2(f,a,b,n)
h = (b-a)/n;
w = h*[0.5,ones(1,n-1),0.5];
x = linspace(a,b,n+1);
I = sum(w.*f(x));
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
