diff --git a/int_prob.m b/int_prob.m
new file mode 100644
index 0000000000000000000000000000000000000000..79cf09cb3eaf16094242a826a65c852982364613
--- /dev/null
+++ b/int_prob.m
@@ -0,0 +1,29 @@
+function p = int_prob(l, T, x0, direc, ind_t, bp, ind_delta, T_mov)
+% direc ... change direction of jumps, 1: left->right, -1: right->left
+% ind_t     ... time index at which propagators are evaluated
+% bp        ... boundary position at t==ind_t
+delta_x0 = diff(x0);
+p = 0;
+for i = 1:length(delta_x0)
+    x = (x0(i)+x0(i+1))/2;
+    % 1. cond.: corr. starting point?  2. cond: jumped outside of domain?
+    corr_starting_point = direc*x < direc*bp;
+    corr_end_point = direc*(x-l) > direc*(bp-T{i}.v*T{1}.t(ind_delta+1));
+    if corr_starting_point && corr_end_point
+        if nargin==7
+            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_t, :), x-l);
+            p2 = (p_i(i)+p_i(i+1))/2;
+            % @ SS distribution for x0 can be taken from phi_tot
+            p = p + delta_x0(i)*...
+                    T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0)*p2;
+        elseif nargin==8
+            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :), x-l);
+            p2 = (p_i(i)+p_i(i+1))/2;
+            p_sol = @(j) interp1(T_mov.x, T_mov.sol(ind_t, :), x);
+%             p_sol2 = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0);
+            p_sol2 = (p_sol(i)+p_sol(i+1))/2;
+            p = p + delta_x0(i)*p_sol2*p2;
+        end
+    end
+end
+end
\ No newline at end of file
diff --git a/prob_laplace.m b/prob_laplace.m
index 58e27746b6aad5f1410f5928b2c45fdf5f8c546e..7c814d56f4ac532e784816baa4b19513c9187580 100644
--- a/prob_laplace.m
+++ b/prob_laplace.m
@@ -417,35 +417,7 @@ csvwrite('prob_in_out.csv', p);
 
 %% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function p = int_prob(l, T, x0, direc, ind_t, bp, ind_delta, T_mov)
-% direc ... change direction of jumps, 1: left->right, -1: right->left
-% ind_t     ... time index at which propagators are evaluated
-% bp        ... boundary position at t==ind_t
-delta_x0 = diff(x0);
-p = 0;
-for i = 1:length(delta_x0)
-    x = (x0(i)+x0(i+1))/2;
-    % 1. cond.: corr. starting point?  2. cond: jumped outside of domain?
-    corr_starting_point = direc*x < direc*bp;
-    corr_end_point = direc*(x-l) > direc*(bp-T{i}.v*T{1}.t(ind_delta+1));
-    if corr_starting_point && corr_end_point
-        if nargin==7
-            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_t, :), x-l);
-            p2 = (p_i(i)+p_i(i+1))/2;
-            % @ SS distribution for x0 can be taken from phi_tot
-            p = p + delta_x0(i)*...
-                    T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0)*p2;
-        elseif nargin==8
-            p_i = @(j) interp1(T{j}.x, T{j}.sol(ind_delta+1, :), x-l);
-            p2 = (p_i(i)+p_i(i+1))/2;
-            p_sol = @(j) interp1(T_mov.x, T_mov.sol(ind_t, :), x);
-%             p_sol2 = T{1}.phi_tot(x, T{1}.a, T{1}.b, T{1}.e, T{1}.u0, 0);
-            p_sol2 = (p_sol(i)+p_sol(i+1))/2;
-            p = p + delta_x0(i)*p_sol2*p2;
-        end
-    end
-end
-end
+
 
 function p = normalization(T, x0, lb, ind_t, direc, bp)
 delta_x0 = diff(x0);