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 )
y1=ic(1);
y2=ic(2);

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

Calculating Solutions:

clear; close all; clc
set(0,'defaultlinelinewidth',2);
set(0,'defaultaxeslinewidth',1.5);
set(0,'defaultpatchlinewidth',2);
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
dbeta=1;
betasol=[];
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;
else
beta = beta + dbeta/2; %uses bisection to converge to the solution
dbeta=dbeta/2;
end
if abs(y(end,1)) < tol % check for convergence
betasol=[betasol beta]; % write out eigenvalue
break % get out of convergence loop
end
end
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
end
betasol
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')
saveas(gcf,'DiffEqnSol','jpg')