function varargout=shoot(odefun,xspan,A,B,varargin)
%SHOOT Solves boundary-valued second-order ODEs using a shooting method.
%   [X,U] = SHOOT(G,XSPAN,A,B) with XSPAN = [X0, X1, X2, ..., XN] is the 
%   solution of the second-order ODE U''(X) = F(X,U(X),U'(X)), with boundary
%   values U(X0) = A, U(XN) = B where G is a function handle for the associated
%   system of first-order ODEs. The equation is integrated by stepping from 
%   X0 to X1 to XN using ODE4; if XSPAN = [X0 XN] then ODE45 is used with 
%   RelTol = 1e-6 and AbsTol = 1e-6.
%
%   The method consists in solving iteratively the initial value problem 
%   V''(X) = F(X,V(X),V'(X)), V(X0) = P, V'(X0) = Y, where Y is an intially 
%   guessed slope at X = X0, until the root of the function E(Y) = V(X0) - B 
%   is found using a secant method.
%
%   [X,U] = SHOOT(G,XSPAN,A,B,P1,P2,...) passes the parameters P1,P2,... to G.
%
%   [X,U,Y,E] = SHOOT(...) returns E(Y).
%
%   See also SECANT.

%   Francisco J. Beron-Vera, 2006/11/05

if length(xspan)==2
   x=linspace(xspan(1),xspan(2),25);
   options = odeset('RelTol',1e-6,'AbsTol',1e-6);
else
   x=xspan;
end
N = 25;
u = zeros(length(x),N+2);
y = zeros(N+2,1);
E = zeros(N+2,1);

y(1) = (B-A)/(xspan(end)-xspan(1));
if length(xspan) == 2
   [x z] = ode45(odefun,x,[A;y(1)],options,varargin{:});
else
   z = ode4(odefun,x,[A;y(1)],varargin{:});
end
u(:,1) = z(:,1);
E(1) = z(end,1)-B;

if sign(E(1)) > 0
   y(2) = .5*y(1);
else
   y(2) = 2*y(1);
end
if length(xspan) == 2
   [x z] = ode45(odefun,x,[A;y(2)],options,varargin{:});
else
   z = ode4(odefun,xspan,[A;y(2)],varargin{:});
end
u(:,2) = z(:,1);
E(2) = z(end,1)-B;

n = 1;
while abs(y(n+1)-y(n)) > eps*abs(y(n+1))
   n = n+1;
   if n > N
      warning('Tolerance not satisfied.')
      break
   end
   y(n+1) = y(n)-E(n)*(y(n)-y(n-1))/(E(n)-E(n-1));
   if length(xspan) == 2
      [x z] = ode45(odefun,x,[A;y(n+1)],options,varargin{:});
   else
      z = ode4(odefun,x,[A;y(n+1)],varargin{:});
   end
   u(:,n+1) = z(:,1);
   E(n+1) = z(end,1)-B;
end

if length(xspan)==2
   varargout{1} = x;
else
   varargout{1} = x;
end
varargout{2} = u(:,1:n);
if nargout > 2
   varargout{3} = y(1:n);
   varargout{4} = E(1:n);
end