function Shooting



%%  Use MATLAB's nonlinear equation solver to compute parameter Alpha for
%%  which the 1st component y1 of the solution to the ODE IVP
%%    y1' = y2
%%    y2' = r(t), 0 < t < 1
%%    y1(0) = 0
%%    y2(0) = Alpha
%%  matches the solution to the ODE BVP
%%      u'' = r(x), 0 < x < 1,
%%      u(0) = 0,  u(1) = 0.
%%  Note that alpha is a MATLAB function; hence the choice of Alpha for the
%%  left slope parameter.

%%  Set solver options.
fsolve_options = optimset('Display','iter','LargeScale','off');

%%  Set initial guess for slope u'(0) = Alpha.
Alpha0 = 0;

%%  Solve equation to obtain Alpha for which y(1;Alpha) matches u(1) = 0.
Alpha = fsolve(@ODEeval,Alpha0,fsolve_options);

%%  Display results.
figure(1)
  plot(tvec,yvec(:,1))
  xlabel('x')
  ylabel('u(x)')
  title('solution to ODE BVP')

%%  Nested function solves ODE initial value problem with specified slope
%%  at t = 0 and returns solution value at t = 1.

    function u_R = ODEeval(slope)
        
        yvec_init = [0; slope];
        time_interval = [0,1];
        
        %%  Compute (constant) Jacobian.
        J = [0 1; 0 0];
  
        %%  Set ODE solver options.
        ODEoptions = odeset('Jacobian',J);
       
        %%  Solve ODE IVP.
        [tvec,yvec] =  ode15s( @rhs, time_interval, yvec_init, ODEoptions );
        
        %%  Extract value of solution at t = 1.
        u_R = yvec(end,1);
        
    end

    function fvec = rhs(t,y)
        
        fvec = zeros(2,1);
        fvec(1) = y(2);
        fvec(2) = -1;  % take r(t) = -1
        
    end

end