function [tout, yout] = thetametodo (odefun, tspan, y0, options)
%
% [tout, yout] = thetametodo (odefun, tspan, y0, options)
%
% Solve the ODE y'(t) = odefun (t, y(t)), where
%
% odefun is ...
% tspan is ...
% y0 is ...
% tout is ...
% options.Theta is ...
m = length (tspan) - 1;
yout = zeros (length (y0), m + 1);
yout(:, 1) = y0(:); % make y0 a column vector
k = tspan(2) - tspan(1);
theta = options.Theta;
dodefundy = options.Jacobian;
if (theta == 0)
  error ('Do not use thetametodo with theta = 0')
end
tol = k^2 / 100;
for n = 1:m
  F = @(x) x - yout(:,n) - k * (1 - theta) * odefun (tspan(n), yout(:,n)) - ...
           k * theta * odefun (tspan(n + 1), x);
  J = @(x) eye (length (x)) - k * theta * dodefundy (tspan (n+1), x);
  yout(:, n+1) = yout(:, n); % initial guess
  delta = -J(yout(:,n+1)) \ F(yout(:,n+1));
  while (norm (delta) > tol) % absolute error estimate
    yout(:,n+1) = yout(:,n+1) + delta;
    delta = -J(yout(:,n+1)) \ F(yout(:,n+1));
   end
  yout(:,n+1) = yout(:,n+1) + delta;
end
tout = tspan(:);
yout = yout.';

