function [t,x]=traprule(f,x0,ti,tf,h)
  t=[ti:h:tf];
  x=zeros(length(x0),length(t));
  x(:,1)=x0;
  for k=1:length(t)-1
    F=@(z) z-(x(:,k)+h/2*(f(z,t(k+1))+f(x(:,k),t(k))));
    x(:,k+1)=fsolve(F,x(:,k));
  end
end