function [oneshot, shooting] = bvpshoot(drag,freq,forcing,ybcleft,ybcright)
%
%  bvp_shoot returns two function handles that help solve the bvp:
%            y''+drag(t)*y'+freq(t)*y = forcing(t),   
%                 y(0)=ybcleft,  y(10)=ybcright
%
%    inputs:      
%       drag,freq,forcing are function handles for the coefficients
%       ybcleft and ybcright are the scalar constant y-values at the boundary.
%
%    outputs:
%       oneshot--function handle for shooting across once
%             inputs:  guess--the velocity at the left boundary
%             outputs: result--the difference between y(10) and ybcright
%                      t,y--the vector of times and y-values for the shot
%       shooting--function handle for shooting once a bracket is found.
%       
%       example usage: solve y''+2*y'+3*y=1, y(0)=2,  y(10)=4 
%          dr=@(t) 2;
%          fr=@(t) 3;
%          fo=@(t) 1;
%          [shot, shooting]=bvp_shoot(dr,fr,fo,2,4);
%          result1=shot(0);
%          result2=shot(1); % Use "shot" to bracket the root.
%          [t,y]=shooting(0,1)  % Give the bracket to the shoting routine.
%          plot(t,y(:,1))
%
  function [result,t,y]=shot(guess)
    yinit=[ybcleft, guess];
    rhs=@(t,y) [y(2);-drag(t).*y(2)-freq(t).*y(1)+forcing(t)];
    tspan=[0 10];
    [t,y]=ode45(rhs,tspan,yinit);
    result=y(end,1)-ybcright;
  end
  oneshot=@shot;

  function [t, y]=shoot(a,b)
    root=bisection(oneshot,a,b);
    [r,t,y]=oneshot(root);
  end
  shooting=@shoot;
end
