clear all;
close all;

% Parameters of the system:
a = 1;
b = 5;
h = input('Enter the stepsize h = ');
r = a : h : b;
n = length(r)-1;

% Matrix set-up for numerical solution based on central differences:
M = 2*diag(ones(1,n-1))-diag(ones(1,n-2),1)-diag(ones(1,n-2),-1) ...
    -0.5*h*diag(1./r(2:n-1),1) + 0.5*h*diag(1./r(3:n),-1);
[V,D] = eig(M);
[omega,index] = sort(sqrt(diag(D)/h^2));
smallesteigevalues = omega(1:5)

% Plotting eigenvectors for the first smallest eigenvalues
Vector = V(:,index(1:5)); Vector = [ zeros(1,5); Vector; zeros(1,5) ]; % adding boundary values to eigenvectors

figure(1); 
plot(r',Vector(:,1)); xlabel('r'); ylabel('first eigenvector'); 
figure(2); 
plot(r',Vector(:,2)); xlabel('r'); ylabel('second eigenvector'); 
figure(3); 
plot(r',Vector(:,3)); xlabel('r'); ylabel('third eigenvector'); 
figure(4); 
plot(r',Vector(:,4)); xlabel('r'); ylabel('fourth eigenvector'); 
figure(5); 
plot(r',Vector(:,5)); xlabel('r'); ylabel('fifth eigenvector'); 
