clear all
close all
%setting variables
walk_num = 100;
W = 0.1;
N = 100000;
beta = 1;

%assingment 12: plotting the potencial
x = linspace(-3,3,100);
figure()
plot(x,V(x),'o-')
set(gca,"fontname","arial","fontsize",14,"fontweight","normal");
ylabel ("\it{V(x)}");
xlabel ("\it{x}");
grid on
print("figure_4_1", "-dpng", "-r300");


%assingment 14
new_walkers = zeros(walk_num,1);
old_walkers = -3 + 6*rand(walk_num,1); %initalising walkers
for n = 1:N %monte carlo simulation
p_old = exp(-beta.*V(old_walkers));
new_walkers = old_walkers+(-W + 2.*W.*rand(walk_num,1));
p_new = exp(-beta.*V(new_walkers));
ratio = p_new./p_old - rand();
new_walkers(ratio<0) = old_walkers(ratio<0);
old_walkers = new_walkers;
end


figure()
counts = histcounts(new_walkers,x);
histogram('BinEdges',x,'BinCounts',counts/walk_num) %plotting numerically values
hold on
p = exp(-beta.*V(x))/(sum(exp(-beta.*V(x))));
plot(x,p) %potting theoretical values



grid on
set(gca,"fontname","arial","fontsize",14,"fontweight","normal");
xlabel ("\it{x}");
ylabel ("Probability");
print("figure_4_2", "-dpng", "-r300");



function [v] = V(x)
v=x.^2+2.*exp(-4.*(x-0.5).^2)+2.5.*exp(-4.*(x+0.5).^2);
end