% Shooting method with bisection
% 5-20-2008
%
% Note that you will need the function bvp_shooting_error.m and 
% bvp_shooting_rhs.m in order for this script to run!
%
% Shooting method for a simple boundary value problem
% Problem statement:
% y''(t) = -sin(y(t))   (Nonlinear pendulum)
% y(0) = 0, y(2) = pi/2
% Find y(t)

clear all; close all; clc

TOLERANCE = 10^(-4);  % Tolerance for the bisection
max_iterations = 100;  % Maximum number of bisection steps
plot_iterations = false;  % Plot each iteration for debugging

tspan = [0 2];  % Domain

A = 1;  % Initial guess for the slope at the rhs
dA = 0.5;  % How much to initialy change A

for i=1:max_iterations
    y0 = [0 A];  % Construct the initial condition
    [t,y] = ode45('bvp_shooting_rhs',tspan,y0);  % Solve the ODE
    
    if abs(y(end,1) - pi/2) < TOLERANCE  % Check to see if we are done
        break;
    end
    if y(end,1) < pi/2  % Modify A according to where we are wrt pi/2
        A = A + dA;
    else
        A = A - dA;
        dA = dA/2;  % Half our step
    end
    
    if plot_iterations  % Plot each iteration if plot_iterations = true
        hold on;
        plot(t,y(:,1));
        plot(2,pi/2,'r.','MarkerSize',15);
        hold off;
    end
end
    
plot(t,y(:,1),'rx'); hold on;

% Compare to our fminsearch approach
yp0 = fminsearch('bvp_shooting_error', 1);

y0 = 0;
[t,y] = ode45('bvp_shooting_rhs', [0 2], [y0 yp0]);

plot(t, y(:,1), 'b-')

title('Solution Comparison'); xlabel('x'); ylabel('y(x)');
legend('Bisection','fminsearch','Location','SouthEast');