function [tvec,yvec] = Trapezoid(t0,y0,f,h,N,tol)
% [tvec,yvec] = Trapezoid(t0,y0,f,h,N)
% Trapezoidal rule
% Inputs
% t0,y0: initial condition
% f:     names of the right-hand side function f(t,y)
% h:     stepsize
% N:     number of steps
% tol:   tolerance for fixed point iteration
% Outputs
% tvec: vector of t values
% yvec: vector of corresponding y values

yvec = zeros(N+1,1);
tvec = linspace(t0,t0+N*h,N+1)';
yvec(1) = y0;
w = yvec(1);
for n = 1:N
   fvalue1 = f(tvec(n),w);
   converged = 1;   k = 1;
   while  (converged > tol) & (k < 10)     % fixed point iteration to solve implicit problem
       fvalue2 = f(tvec(n+1),w);
       w = yvec(n)+h/2*(fvalue1+fvalue2);
       if k > 1
           converged = abs(wold-w);
       end
       wold = w;
       k = k+1;
   end    
   yvec(n+1) = w;
end