clc,clear all,close all
dx=1;L=500;tmax=400;
dt=2;x=0:dx:L;I=length(x);
%% condição inicial
sigma=12;x0=180;k0=1;
E=k0^2/2;V0=0.9*E;
mu=1/sqrt(sigma*sqrt(pi));
f0=mu*exp(-(x-x0).^2/(2*sigma^2)).*exp(-i*k0*(x-x0));
f0(1)=0;f0(I)=0;f02=conj(f0').*f0';
%% parametros
hbar=1;m=1;
f=f0';J=tmax/dt;N=I;
%% Matriz do hamiltoniano
D2b=diag(ones([N-1,1]),-1)-2*diag(ones([N,1]),0)+diag(ones([N-1],1),1);
D2b(1,1)=0;D2b(1,3)=0;;D2b(N,N-2)=0;D2b(N,N)=0;
D2b=D2b/dx^2;
%% potencial degrau
V=V0./(1+exp(-150*(x-250)));
%% potencial barreira
%V=V0./(1+exp(-150*(x-250)))-V0./(1+exp(-150*(x-(18+250))));
%% potencial barreira dupla
%V=V0./(1+exp(-150*(x-250)))-V0./(1+exp(-150*(x-(20+250))))+...
%    V0./(1+exp(-150*(x-(250+60))))-V0./(1+exp(-150*(x-(80+250))));
%% potencial harmonico
%V=0.05*(x-L/2).^2;
H=(-hbar^2/(2*m))*D2b+diag(V);
%% Propagando no tempo
t=0;
for j=1:J-1
t=t+dt;
U=expm(-i*H*dt/hbar);
f=U*f;f(1)=0;f(I)=0;
f2=conj(f).*f;

yyaxis left
plot(x,f2,'k-')
xlabel('x','FontSize',18)
ylabel('|\Psi(x,t)|^2','FontSize',18,'Color','k')
set(gca,'FontSize',18)
axis([0 L 0 0.1])
yyaxis right
plot(x,V,'-k');
ylabel('E_p(x)','FontSize',18,'Color','k')
axis([0 L 0 0.6])
set(gca,'FontSize',18)
ax = gca;
ax.YAxis(1).Color = 'k';
ax.YAxis(2).Color = 'k';
pause(.001)
end