% Usage: 
% y = rkf45bisectshoot(f,a,b,ya,yb,s1,s2,tol) or
% [y t] = rkf45bisectshoot(f,a,b,ya,yb,s1,s2,tol) or 
% [y t s] = rkf45bisectshoot(f,a,b,ya,yb,s1,s2,tol)
% Shooting method for boundary value problems using 
% Runge-Kutta-Fehlberg order 4 / order 4 embedded pair.
% Bisection method is employed to solve F(s) = 0
%
% Input:
% f - Matlab inline functions f(t,y,z) where z=y'
% a,b - interval
% ya, yb - boundary conditions
% [s1 s2] - interval containing a zero of F
% tol - tolerance
%
% Output:
% y - computed solution
% t - time steps
% s - value of y'(a)
%
% Examples:
% [y t]=rkf45bisectshoot(inline('-y','t','y','z'),0,1,1,1,0,1,1e-6);

function [y t s] = rkf45bisectshoot(f,a,b,ya,yb,s1,s2,tol)
% transform to a firts order equation in two unknowns
ff = @(t,u) [u(2) f(t,u(1),u(2))];
% initial value of step-size
h = (b-a) / 10;

% compute F(s1)
ua = [ya s1];
[u t] = rkf45(ff,a,b,ua,h,tol);
F1 = u(n+1,1) - yb;
% compute F(s2)
ua = [ya s2];
[u t] = rkf45(ff,a,b,ua,h,tol);
F2 = u(n+1,1) - yb;

if sign(F1)*sign(F2) >= 0
  error('F(s1)F(s2)<0 not satisfied!') %ceases execution
end

while (s2 - s1) / 2 > tol
  sn = (s1 + s2) / 2;
  % compute F(sn)
  ua = [ya sn];
  [u t] = rkf45(ff,a,b,ua,h,tol);
  Fn = u(n+1,1) - yb;

  if Fn == 0              %sn is a solution, done
    break
  end
  if sign(Fn)*sign(F1)<0  %s1 and sn make the new interval
    s2 = sn; F2 = Fn;
  else                    %sn and s2 make the new interval
    s1 = sn; F1 = Fn;
  end
end
% compute output
s = (s1 + s2) / 2;               %new midpoint is best estimate
ua = [ya s];
[u t] = rkf45(ff,a,b,ua,h,tol);
y = u(:,1);