% Script to solve the pendulum ODE BVP
%
%    u" + sin(u) = 0, u(0) = a, u(T) = b
%
% for specified a, b, T using the shooting method. The script is set up to 
% duplicate the parameters chosen for fig. 2.4 of the notes. 
% This implementation uses a  

  a = 0.7;
  b = 0.7;
  T = 2*pi;

%  Make two different starting guesses v0_old and v0_new for the
%  initial value of u'(0), and find u(T) for each of them using the ODE
%  solver ode45.  Then use the secant method to successively update u'(0)
%  until a convergence threshold is reached.  A maximum of 20 iterations is 
%  allowed to prevent an infinite loop.

  v0_old = 0.5;
  v0_new = 0;
  y0_old = [a v0_old];
  tspan = [0 T];
  [t, y_old] = ode45('pendulum_dydt',tspan,y0_old);
  uT_old = y_old(length(t),1);

  subplot(2,2,1)
  plot(t,y_old(:,1),'--')
  axis([0 2*pi -1 1])
  xlabel('t')
  ylabel('\theta (radians)')
  title('Pendulum ODE BVP Shooting Method')
  hold on

  for iter = 1:20
    y0_new = [a v0_new];
    [t, y_new] = ode45('pendulum_dydt',tspan,y0_new);
    uT_new = y_new(length(t),1);
    plot(t,y_new(:,1),'--')

%    Secant method for new guess

    
%    Let upr_new and upr_old be the most and next-most recent guesses for the
%    initial value of u'(0).

    xi =  (uT_new-b)/(uT_new-uT_old);
    delta = xi*(v0_old - v0_new)

%    Exit loop if |delta| < 1e-8 and plot solution

    if(abs(delta) < 1e-8)
      plot(t,y_new(:,1),'-')
      hold off
      v0 = v0_new
      break
    end

%    Update the 'new' and 'old' values of v0

    v0_next = v0_new + delta;
    v0_old = v0_new;
    v0_new = v0_next;
    uT_old = uT_new;

  end
