diff --git a/prob_laplace.m b/prob_laplace.m
index 3e5cb4df635261a74ad3ed4104bdc95eaa1b3110..737bb6f3a7c7b4b27ca3fe19560a4d9161afa8f7 100644
--- a/prob_laplace.m
+++ b/prob_laplace.m
@@ -215,7 +215,11 @@ parfor i = 1:length(ls)
 end
 toc
 %% Normalization factor
-N = normalization(T, x0, 0, 3);
+N = normalization(T, x0, 0, 3, direc, 5)
+% N = sum(p)*0.001;
+% m = sum(ls.*p/N)/3001*3;
+%%
+sum((T{1}.sol(1, 1:end-1)+T{1}.sol(1, 2:end))/2.*diff(T{1}.x))
 %% Do we need to look at the left side as well?
 figure; hold on;
 plot(ls, p/N);
@@ -421,16 +425,16 @@ for i = 1:length(delta_x0)
 end
 end
 
-function p = normalization(T, x0, lb, ind_t)
+function p = normalization(T, x0, lb, ind_t, direc, bp)
 delta_x0 = diff(x0);
-logi = T{1}.x<5-lb;
+logi = -direc*(T{1}.x) < -direc*(bp-lb);
 delta_x = diff(T{1}.x(logi));
 p_i = zeros(1, length(delta_x0));
 for i = 1:length(delta_x0)
     sol = T{i}.sol(ind_t, logi);
     sol = (sol(1:end-1)+sol(2:end))/2;
     x = (x0(i)+x0(i+1))/2;
-    if x > 5+lb
+    if -direc*x > -direc*(bp+lb)
         p_x0i = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0);
         p_i(i) = delta_x0(i)*sum(sol.*delta_x)*p_x0i;
     end