clear; close all; clc;

%%
tol = 1e-6; % Tolerance for the shooting method
xspan = [-1,1];

% Differential equation: d^2y / dt^2 = (beta-100)*y
% Subject to boundary conditions y = 0 at x = +/-1
ic = [0 1]; % Initial conditions
beta = 97;

figure(1); clf;
axes; hold on;
plot(xspan,[0 0],'k--');
xlabel('x'); ylabel('Psi');

for jj = 1:5
    % Iterate through values of beta until the solution satisfies the
    % boundary conditions
    dbeta = 1; % Delta beta
    for j = 1:1000
        odefun = @(t,y) [y(2) (beta-100)*y(1)]';
        [t,y] = ode23(odefun,xspan,ic);
        if abs(y(end,1))<tol
            break;
        end
        if y(end,1)*((-1)^(jj+1)) > 0 % Switch sign of the solution
            beta = beta - dbeta;
        else
            dbeta = dbeta/2;
            beta = beta + dbeta;
        end
        
    end
    fprintf([int2str(jj),'->; beta-100 = ',num2str(beta-100,'%0.5f'),';\n']);
    beta = beta-0.1;
    
    psi = y(:,1);
    norm = trapz(t,psi.^2);
    psi = psi./sqrt(norm);
    
    plot(t,psi,'.-');
    grid on; box on; axis square;
    ylim([-1 1]*1.5);
    drawnow;
end