function  v=odetrap(f,t,dt,u,varargin)
% v=odetrap(f,t,dt,u,varargin)
% one step of the trapezoidal method for u=f(t,u)
%
% arguments
%   f   the rhs function
%   t   starting time
%   dt  time step
%   u   solution at time t
%   varargin    extra arguments to be passed to f
% output
%   v   solution at time t+dt

% update formula:
% v = u + 0.5*dt*(f(t+dt,v)+f(t,u))

% predictor:
v0 = u + dt*feval(f,t,u,varargin{:});
% solve the trapezoidal equation res(du,...)=0
du = nsolve(@res,zeros(size(v0)),v0-u,f,u,t,dt,varargin{:});
% update the solution
v  = u + du;

function r=res(du,f,u,t,dt,varargin)
% residual in the nonlinear system solved in the trapezoidal method
% to solve for increment du the update formula
% v = u + 0.5*dt*(f(t(k+1),v)+f(t,u))
% where v = u + du
f1 = feval(f,t+dt,u+du,varargin{:});
f2 = feval(f,t,u,varargin{:});
r  = du - 0.5*dt*(f1+f2);