function y = euler_midstart(ff,x,x0,y0);
% EULERTAYLOR   Return Euler's method solution for the problem
%       y' = f(x,y),   y(x0) = y0,
%   where x0 is somewhere in the middle of the given x grid.
%   ff.m must compute right-hand side.  Example:
%     >> x = -1:.01:4;
%     >> y = euler_midstart('ff',x,0,2);

hmin = min(diff(x));
i0 = min(find(x>x0-(hmin/2)));

y = zeros(size(x));
y(i0) = y0;
 
%forward from i0
for j=i0+1:length(x)
  y(j) = y(j-1) + (x(j)-x(j-1)) * feval(ff,x(j-1),y(j-1));;
end

%backward from i0
for j=i0-1:-1:1
  y(j) = y(j+1) - (x(j+1)-x(j)) * feval(ff,x(j+1),y(j+1));;
end