function [x, u, dudx] = shoot(phi, xspan, lval, lder, rval, rder, init)
% SHOOT   Shooting method for a two-point boundary-value problem.
%
% Input:
%   phi      defines u'' = phi(x,u,u') (function)
%   xspan    endpoints of the domain (vector)
%   lval     prescribed value for u(a) (use [] if unknown) 
%   lder     prescribed value for u'(a) (use [] if unknown) 
%   rval     prescribed value for u(b) (use [] if unknown) 
%   rder     prescribed value for u'(b) (use [] if unknown) 
%   init     initial guess for lval or lder (scalar)
%
% Output:
%   x        nodes in x (length n+1)
%   u        values of u(x)  (length n+1)
%   dudx     values of u'(x) (length n+1)

    v = [0, 0];  % fake initialization
    % Tolerances for IVP solver and rootfinder.
    ivp_opt = odeset('reltol', 1e-6, 'abstol', 1e-6);
    optim_opt = optimset('tolx', 1e-5);

    % Find the unknown quantity at x=a by rootfinding. 
    s = fzero(@objective, init, optim_opt); 

    % Don't need to solve the IVP again. It was done within the
    % objective function already.
    u = v(:, 1);            % solution     
    dudx = v(:, 2);         % derivative

      % Difference between computed and target values at x=b. 
      function F = objective(s)
          if isempty(lder)   % u(a) given
            v_init = [ lval; s ];
          else               % u'(a) given
            v_init = [ s; lder ];
          end
      
          [x, v] = ode45(@shootivp, xspan, v_init, ivp_opt);
      
          if isempty(rder)   % u(b) given
            F = v(end,1) - rval;
          else               % u'(b) given
            F = v(end,2) - rder;
          end
      end

      % ODE posed as a first-order equation in 2 variables. 
      function f = shootivp(x, v)
          f = [ v(2); phi(x,v(1),v(2)) ];
      end

end