%amonia.m
clear all
dx=0.1; L=5;x=-L:dx:L; N=length(x);
D2=diag(ones([N-1,1]),-1)-2*diag(ones([N,1]),0)+diag(ones([N-1],1),1);
D2=D2/dx^2;
hbar=1;m=1;k=1;b=4;
pot=0.5*k*x.^2+b*exp(-x.^2);H=(-hbar^2/(2*m))*D2+diag(pot);
[V,E]=eig(H);Ediag=diag(E);energ=Ediag(1:8)

figure(1)
plot(x,pot),hold on, 
plot(x,energ(1)),plot(x,energ(2))
plot(x,energ(3)),plot(x,energ(4))
plot(x,energ(5)),plot(x,energ(6))
plot(x,energ(7)),plot(x,energ(8))
plot(x,V(:,1)+energ(1)),plot(x,V(:,2)+energ(2))
plot(x,V(:,3)+energ(3)),plot(x,V(:,4)+energ(4))
plot(x,V(:,5)+energ(5)),plot(x,V(:,6)+energ(6))
plot(x,V(:,7)+energ(7)),plot(x,V(:,8)+energ(8))
xlabel('x'),ylabel('E_p(x),  E,  \psi(x)')
axis([-5,5,min(pot),10])



function f=AjNH3ohgauss(p)
dx=0.001;L=2.5;x=-L:dx:L; N=length(x);n=length(x);EP=p(1)*x.^2+p(2)*exp(-p(3)*x.^2);
D2=toeplitz([-2 1 zeros(1,n-2)])/dx^2;H=-1/(2*4668)*D2+diag(EP);
[V,E]=eig(H);en=diag(E)'; norma=sqrt(trapz(x',V.^2));V=V./norma;
Energia=(en(1:5)-en(1));data=[0 0.7935 932.51 968.32 1597.6]*4.55633528121e-06;
f=norm(data-Energia);


p=[.07598 .05684 1.3696];
x = fminsearchbnd('AjNH3ohgauss',p,[0 0 0],[inf inf inf]);