function [t,y]=mnhf_modeuler(Fun,y0,h,tmax)
%MNHF_MODEULER solves an autonomous ODE using the Modified Euler scheme.
%
% Fun  - the name of the external function
% y0 - initial condition i.e. y(t_0)
% h - step size or time increment
% tmax - maximum time to which to integrate
%
% Usage  mnhf_modeuler(@odefun1,-0.5,1e-2,10.0);

tol = 1e-5; % Set tolerance.
Nmax = 1e3; % Set maximum number of iterations.

% Initialize arrays. Specify initial condition.
t = 0:h:tmax;
y = zeros(size(t)); y(1) = y0;

% Determine number of time steps, N.
N = length(t)-1;

for ij = 1:N
 count = 0;
 y(ij+1) = y(ij)+h*feval(Fun,y(ij));% Explicit Euler provides initial guess
 while abs( y(ij+1)-y(ij)-0.5*h*(feval(Fun,y(ij))+feval(Fun,y(ij+1))) )>tol
  y(ij+1) = y(ij)+0.5*h*(feval(Fun,y(ij))+feval(Fun,y(ij+1)));
  count = count+1;
  if count > Nmax
   error('Maximum number of iterations exceded.\n\n')
  end
 end
end