% quantum mc for harmonic oscillator

rngseed = 32;
rng(rngseed);

n0 = 500000;
x = zeros(n0, 1);

er = 0;
deltat = 0.005;
step = sqrt(deltat);
nsteps = 1000;
nrepls = zeros(nsteps, 1);
generg = zeros(nsteps, 1);

for i=1:nsteps
    lx = length(x);
    nrepls(i) = lx;
    x = x + step*randn(lx, 1); % diffusion step
    v = harmonic(x);
    er = mean(v) - 0.1*(lx/n0 - 1); % added empirical 'feedback' term 
                                    % to keep lx \approx n0 
    generg(i) = er;
    fprintf('%10i  %10i  %.5f\n', i, lx, er);
    w = exp((-v + er)*deltat);
    u = rand(lx, 1);
    m = floor(w + u);
    x2 = x(m >= 2); 
    x3 = x(m >= 3);
    x4 = x(m >= 4);
    x = x(m ~= 0); % death of replicas
    x = [x' x2' x3' x4']'; % birth of replicas
end

[y, edges] = histcounts(x);
ne = length(edges);
mids = (edges(1:ne-1) + edges(2:ne))/2;
boxw = edges(2:ne) - edges(1:ne-1);
y = (y ./ boxw) / sqrt(sum((y .^ 2) ./ boxw)); % normalized wave function
z = 1/sqrt(sqrt(pi)) * exp(-(mids.^2)/2);  % exact wave function

plot(mids, y)
hold on
plot(mids, z)
hold off

function v = harmonic(x)
% Potential for harmonic oscillator
    v = x.^2/2;
end
