We are going to solve the differential equation with the boundary conditions

\psi_{xx} +[100 - \beta ]\psi = 0

\implies \psi_{xx} = [\beta - 100] \psi

Let’s take the simpler boundary condition

\psi(x=\pm 1) =0.

Now, the first thing which we wanna do is to convert this equation in the first order form. After that, we can simply use the 4th order Runge-Kutta to obtain higher order solution.

To transform our 2nd order differential equation into first order, we define two variables. To transform nth order equations we need n such variables. This is the standard process to convert any differential equation into the system of first order differential equations.

Lets take

y_1 = \psi (1)

\implies y_1' = \psi_x

\implies y_1' = y_2

y_2= \psi_x (2) [derivative of \psi]

\implies  y_2' = \psi_{xx}

\implies  y_2' = [\beta - 100] \psi

So, we have transformed the 2nd order differential equation into two first order differential equations.

And our boundary condition will become y_1(\pm 1) =0

Now for the differential equations of the type

\psi_{xx} = [\beta - 100] \psi

\implies \psi_{xx} = a \psi

have the solution \psi = e^{\pm \sqrt{a}x} (for positive a). Now this solution doesn’t satisfy our condition \psi(x=\pm 1) =0. Exponential functions are monotonic and hence cannot be 0 at two points. So, this tells us that a is not positive, i.e., \beta \ngtr 100.

Now, if \beta < 100, then our eqn becomes

\psi_{xx} = - a \psi

The solution to this is of the form of cosines and sines, which has potential to satisfy the boundary condition. So, we have \beta < 100.

Defining Function in MATLAB:

function [ rhs ] = yfunc( x,ic,dummy,beta )

rhs = [y2; (beta -100)*y1];


Calculating Solutions:

clear; close all; clc
set(0, 'DefaultAxesFontSize', 14)

tol=10^(-4); % define a tolerance level to be achieved
% by the shooting algorithm

xspan=[-1 1]; %boundary conditions % define the span of the computational domain

A=1; % define the initial slope at x=-1
ic=[0 A]; % initial conditions: x1(-1)=0, x1?(-1)=A

n0=100; % define the parameter n0
beta_start=n0; % beginning value of beta
for jj=1:5 % begin mode loop
    beta=beta_start; % initial value of eigenvalue beta
    dbeta=n0/100; % default step size in beta
    for j=1:1000 % begin convergence loop for beta
        [t,y]=ode45('yfunc',xspan, ic,[],beta); % solve ODEs
       if y(end,1)*((-1)^(jj+1)) > 0 %this checks if the beta needs to be higher or lower for convergence
           beta = beta-dbeta;
           beta = beta + dbeta/2; %uses bisection to converge to the solution
       if abs(y(end,1)) < tol % check for convergence
           betasol=[betasol beta]; % write out eigenvalue
           break % get out of convergence loop
    beta_start=beta-0.1; % after finding eigenvalue, pick
    % new starting value for next mode
    norm=trapz(t,y(:,1).*y(:,1)); % calculate the normalization
    plot(t,y(:,1)/sqrt(norm)); hold on
legend(sprintf('\\beta_1 = %.4f',betasol(1)),sprintf('\\beta_2 = %.4f',betasol(2)),sprintf('\\beta_3 = %.4f',betasol(3)),...
    sprintf('\\beta_4 = %.4f',betasol(4)),sprintf('\\beta_5 = %.4f',betasol(5)));
xlabel('x'); ylabel('\psi (x)','FontWeight','bold')