How to plot the different vibrational wave functions for the ground state of 1D harmonic oscillator?

%Here in this code I want to plot the different vibrational wave functions for the ground electronic state of Iodine molecule.
%I took the potential to be harmonic. I put all the parameters for the iodine molecule from literature. but, unable to plot the different wave functions. 
%for every 'n' in this code, the plot is showing same as the ground vibrational state which is a gaussian. kindly help and suggest!

%The form of the wave function is attached. check the attached SS.   

%%%%%%All the calculations are done in AU (atomic unit)


clc
clear all 

% define various conversion constants
AU_TO_EV=27.2113961;
NM_TO_AU=45.56335;
AU_TO_CM=2.194746e5;
AU_TO_ANG=0.529177249;
AU_TO_PS=0.0000241889;
AMU_TO_AU=1.8228885e3;
conv_r= AU_TO_ANG;
conv_en=AU_TO_CM;
conv_t= AU_TO_PS;


r=(1:0.01:4).'/AU_TO_ANG; %coordinate
dr = r(2) - r(1);      % Coordinate step.

% mass of I2
mass= (126.904473/2)*AMU_TO_AU;
hbar=1;
% X (ground) state of I2---not used here for dynamics but useful for
% visualization 
De_X= 12440/AU_TO_CM;
be_X= 1.875*AU_TO_ANG;
re_X= 2.656/AU_TO_ANG;
we_X= be_X*sqrt(2*De_X/mass);

V1= 1/2*mass*(we_X*(r-re_X).^2);  %ground state potential
% V1= De_X*(1-exp(-be_X*(r-re_X))).^2;

% % initial ground state wavefunction (ground state of harmonic oscillator---an
% % approximation here, but a good one.)
w_x_vib= (212)/AU_TO_CM;
beta= sqrt(mass*w_x_vib);
X=beta*r;
n=2;      %quantum no
Nn=beta./sqrt(sqrt(pi)*(2.^n)*factorial(n));
Nn1=Nn*hermiteH(n,X);
psi_init= Nn1.*exp(-0.5*(beta^2)*((r-re_X).^2));
norm= 1/sqrt(dot(psi_init, psi_init));
Psi_0=psi_init*norm;

plot(r*conv_r,psi_init)