% 5-13-2008
% Direct solve of a BVP

clear all; close all;
n=200;  % Number of grid points
x=linspace(-1,1,n);  % Linearly space out the grid points in [-1,1]
dx=2/n;

n0 = 100;  % Parameter from the ode
m=n-2;  % Size of the matrix we will construct

% Construction of the matrix using for loops
for i=1:m % Main Diagonal
    A1(i,i) = (-2+dx^2*n0)/dx^2;
end
for i=1:m-1 % Off diagonals
    A1(i,i+1) = 1/dx^2;
    A1(i+1,i) = 1/dx^2;
end

% We can also construct the same matrix with one command. This is a
% significantly faster method than using the loops
A2 = 1/dx^2*(diag((-2+dx^2*n0)*ones(1,m),0) + diag(ones(1,m-1),1) ...
    + diag(ones(1,m-1),-1));

% Solve the eigen value problem
[V,D] = eig(A2);

% Plot the first five eigen modes and the eigen values we found
col = ['r','b','k','g','y'];
subplot(1,2,1); hold on;
for mode=1:5
    phi = [0 V(:,end-mode+1)' 0];
    norm = trapz(x,phi.^2);
    phi = phi / sqrt(norm);
    plot(x, phi, col(mode));
end
title('Eigen Functions \Psi_n'); xlabel('x'); ylabel('\Psi_n')
subplot(1,2,2)
plot(diag(D),'o')
title('Eigen Values \beta_n'); xlabel('n'); ylabel('\beta_n')