N = 2.^[1:15];

function S = trapez (fun,a,b,N)
  h = (b-a)/N;
  fy = feval(fun,linspace(a,b,N+1));
  fy(1) = fy(1)/2;
  fy(N+1) = fy(N+1)/2;
  S = h*sum(fy);
end

function y = f(x)
  y = exp(x);
end

function y = g(x)
  y = sin(x*2*pi);
end

function y = h(x)
  y = x.^(1/3);
end

for k=1:15;
  f_err(k) = abs(exp(1)-1-trapez('f',0,1,2^k));
  g_err(k) = abs(trapez('g',0,1,2^k));
  h_err(k) = abs(3/4-trapez('h',0,1,2^k));
end

title('Trapezoidal rule');
xlabel('N');
ylabel('Error');
loglog(N,1./(N.^2),';Slope 2;', ...
       N,f_err,'x;f(x)=exp(x);', ...
       N,g_err,'o;g(x)=sin(2*pi*x);', ...
       N,h_err,'*;h(x)=x^(1/3);')