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 = sinh(x);
end

function y = g(x)
  y = cosh(x);
end

for k=1:15;
  f_err(k) = abs(trapez('f',-1,1,2^k));
  g_err(k) = abs(2*sinh(1)-trapez('g',-1,1,2^k));
end

title('Trapezoidal rule');
xlabel('N');
ylabel('Error');
loglog(N, 1./(N.^2), ';Slope 2;', ...
       N, f_err, 'x;f(x)=sinh(x);', ...
       N, g_err, 'o;g(x)=cosh(x);')