function [t,y]=nlbound(funcn,funcv,tspan,x0f,tol,varargin)
   t0=tspan(1);tfinal=tspan(2);  ga=x0f(1); gb=x0f(2); m=1; m0=0;
   while (norm(m-m0)>tol), m0=m;
      [t,v]=ode45(funcv,tspan,[ga;m;0;1],varargin);
      m=m0-(v(end,1)-gb)/(v(end,3));
   end
   [t,y]=ode45(funcn,tspan,[ga;m],varargin);