function [fun, x] = LinearShooting(pFun, qFun, rFun, BV)
% On input: 
%   
% pFun, qFun, rFun are functions in
% y'' = p(x)*y' + q(x)* y + r(x)
% BV is a data structure containing
% BV.a, BV.b: end points of interval
% BV.alpha, BV.beta: boundary values at end points.
%
% On output:
%   x:  points in the interval.
% fun: function values at points x.
%
% Written by Ming Gu for Math 128B, Spring 2010
% 
[pFunFcn,msg] = fcnchk(pFun,0);
if ~isempty(msg)
    error('InvalidFUN',msg);
end
[qFunFcn,msg] = fcnchk(qFun,0);
if ~isempty(msg)
    error('InvalidFUN',msg);
end
[rFunFcn,msg] = fcnchk(rFun,0);
if ~isempty(msg)
    error('InvalidFUN',msg);
end
y0 = [BV.alpha;0;0;1];
[Tout, Yout] = ode45(@(x,y)ODEFunfromLinearShooting(x,y,pFunFcn, qFunFcn,rFunFcn),[BV.a,BV.b],y0);
x   = Tout(:);
tau = (BV.beta-Yout(end,1))/Yout(end,3);
fun = Yout(:,1)+tau*Yout(:,3);

