diff --git a/prob_laplace.m b/prob_laplace.m
index 549a8fa019f20e07799489a6410661636aa83f62..bc1139a9af60d469eae28c67af339fd66aeacfe0 100644
--- a/prob_laplace.m
+++ b/prob_laplace.m
@@ -402,9 +402,30 @@ end
 params = {-5, b(7/3, 10^-6), 0.5, e(7/3), 0, 1, 10, 7, 0, 'Constituent', 0};
 t = linspace(0, 5, 51);
 bp = 5;
-direc = 1; % 1: jump from left to right, -1: jump from right to left.
+direc = -2; % 1: jump from left to right, -1: jump from right to left.
+
+%% 
+[ls, p_in_in] = frap_distro(t, params, 2, bp);
+[~, p_out_out] = frap_distro(t, params, -2, bp);
+[~, p_in_out] = frap_distro(t, params, 1, bp);
+[~, p_out_in] = frap_distro(t, params, -1, bp);
+%%
+figure; hold on;
+plot(ls, sum(p_in_in, 2));
+plot(ls, sum(p_out_out, 2));
+plot(ls, sum(p_in_out, 2));
+plot(ls, sum(p_out_in, 2));
+%% Save time points for comparison with python
+csvwrite('lb_in_out.csv', ls)
+csvwrite('prob_in_out.csv', p);
+
+%% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [ls, p] = frap_distro(t, params, direc, bp)
 x0 = sort(5-abs(direc)/direc*(0:0.002:4.01));
-%% Run simulations for 'delta' IC across outside
+
+% Run simulations for 'delta' IC across outside
 F = {};
 parfor i = 1:length(x0)
     tic
@@ -413,7 +434,8 @@ parfor i = 1:length(x0)
     F{i}.solve_tern_frap();
     toc
 end
-%%  Calc. probs. for each jump length in ls and sum over time
+
+%  Calc. probs. for each jump length in ls and sum over time
 T_mov = Ternary_model(0, 'FRAP', params, t, 0.2);
 T_mov.solve_tern_frap();
 ls = -abs(direc)/direc*(0.001:0.04:4);
@@ -427,12 +449,7 @@ for j = 1:n_T
     toc
 end
 % save prob_laplace_X_7_3_FRAP_in_out
-%% Save time points for comparison with python
-csvwrite('lb_in_out.csv', ls)
-csvwrite('prob_in_out.csv', p);
-
-%% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+end
 
 
 function p = normalization(T, x0, lb, ind_t, direc, bp)
@@ -452,6 +469,7 @@ end
 p = sum(p_i);
 end
 
+
 function [dudx, x, sol, flux, g0, pt] = calc_flux(Temp, t_ind)
 
 dudx = diff(Temp.sol(t_ind, :))./diff(Temp.x);