function [sVals, xs, eCode] = dirichletShooter(odefn, xspan, ...
			diriBCs, startQ, N, tol, maxIter, linear)
%
%  Inputs:
%         odefn - an inline function of (x, y) that gives slopes y'
%                 of the first-order system y' = f(x, y)
%         xspan - [xfirst xlast]
%       diriBcs - vector [bdryfirst bdrylast] giving specified values
%                 of y at ends of interval
%        startQ - vector [q1 q2] of first two trial values for y'(xfirst)
%             N - number of subintervals into which to break xspan
%           tol - acceptable size for |y(last) - bdrylast|
%       maxIter - largest number of iterations of secant method allowed
%        linear - boolean (0 = false, nonzero = true) indicating if
%                 f is linear in f(x, y)
%
% Outputs:
%            xs - vector of x-values at which final solution has been found
%         sVals - 2-column vector of y, y' values corresponding to xs
%         eCode - reason routine stopped
%                   0 if acceptable soln within tol was reached
%                   1 if maximum number of iterations
%

if (nargin < 8)
  linear = 0;
elseif (linear ~= 0)
  linear = 1;
end

if (linear == 0)
  q1 = startQ(1); q2 = startQ(2);
  bcStart = diriBCs(1); bcEnd = diriBCs(2);
  [sVals, xs] = rk4sys(odefn, xspan, [bcStart; q1], N);
  s = size(sVals);
  e1 = sVals(1, s(2)) - bcEnd;
  eCode = 0;
  q = [q1; q2];

  for j = 0:maxIter
	[sVals, xs] = rk4sys(odefn, xspan, [bcStart; q2], N);
	s = size(sVals);
	e2 = sVals(1, s(2)) - bcEnd;
	if abs(e2) < tol
	  return
	end
	newq = q2 - e2*(q2 - q1)/(e2 - e1);
	q = [q; newq];
	q1 = q2;
	e1 = e2;
	q2 = newq;
  end
  eCode = 1;

else
  q1 = 0; q2 = 1;
  bcStart = diriBCs(1); bcEnd = diriBCs(2);
  [sVals1, xs] = rk4sys(odefn, xspan, [bcStart; q1], N);
  sVals2 = rk4sys(odefn, xspan, [bcStart; q2], N);
  lambda = (bcEnd - sVals2(1,N+1))/(sVals1(1,N+1) - sVals2(1,N+1));
  sVals = lambda*sVals1 + (1 - lambda)*sVals2;
  eCode = 0;
end