function [position,potentialEnergy,wavefunc,eigenEnergy] = ...
obtainVibrationalWavefunction(peinpu,hbar,mass,maxNumber,steplength)
%
position = peinpu(:,1);
potentialEnergy = peinpu(:,2);
potentialOperator = diag(potentialEnergy);
% Discretize kinetic energy
kineticOperator = -2*eye(length(position));
ok1 = eye(length(position));
kineticOperator(1:end-1,:) = kineticOperator(1:end-1,:) +ok1(2:end,:);
kineticOperator(2:end,:) = kineticOperator(2:end,:) +ok1(1:end-1,:);
kineticOperator = -hbar^2/2/mass/steplength^2*kineticOperator;
[wavefunc,eigenEnergy] = eig(potentialOperator+kineticOperator);
eigenEnergy = diag(eigenEnergy(1:maxNumber,1:maxNumber));
wavefunc = wavefunc(:,1:maxNumber);
for ii = 1:maxNumber
wavefunc(:,ii) = wavefunc(:,ii)/sqrt(trapz(position,wavefunc(:,ii).^2));
end
end