% ============================================================================================
% =====           Solving Boundary Value Problems With MATLABs ODE45 Function            =====
% =====                         Using The Shooting Method                                =====
% ===== The Initial Higher Order Equation is:                                            =====
% =====                                         (d^2u(n)/dx^2)+[n(x)-b(n)]u(n)=0         =====
% ============================================================================================

clear
clc

ToleranceLevel=10^(-4);
EigenfunstionColours=['r','b','g','c','m','k'];
Parametern0=100;
InitialSlope=1;
InitialConditions=[0 InitialSlope];
ComputationalDomain=[-1 1];

%Main Program
BetaBeginningValue=Parametern0;

for Modes=1:5
   
   EigenValueBeta=BetaBeginningValue;
   BetaStepSize=Parametern0/100;
   
   for j=1:1000
      
      [TimeInterval,Solution]=ode45('ShootingMethod',ComputationalDomain,InitialConditions,[],Parametern0,EigenValueBeta);
      
      if abs(Solution(end,1)-0)<ToleranceLevel
         
         EigenValueBeta
         break
         
      end
      if (-1)^(Modes+1)*Solution(end,1)>0
         
         EigenValueBeta=EigenValueBeta-BetaStepSize;
         
      else
         
         EigenValueBeta=EigenValueBeta+BetaStepSize/2;
         BetaStepSize=BetaStepSize/2;
         
      end
   end
   
   BetaBeginningValue=EigenValueBeta-0.1;
   Normalization=trapz(TimeInterval,Solution(:,1).*Solution(:,1));
   plot(TimeInterval,Solution(:,1)/sqrt(Normalization),EigenfunstionColours(Modes));
   title('Five Normalized Eigenmodes Along Their Corresponding EigenValues, n0=100');
   xlabel('X');
   ylabel('Y(X)');
   hold on
   grid on
   
end