diff --git a/prob_laplace.m b/prob_laplace.m
index 55325cf5da26283864658e7662df0fb080954455..bb79134e2ecad201745546880b9ca79e3477b160 100644
--- a/prob_laplace.m
+++ b/prob_laplace.m
@@ -348,7 +348,7 @@ end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 t = linspace(0, 5, 300);
 T_mov = Ternary_model(0, 'FRAP', {-5, b(7/3, 10^-6), 0.5, e(7/3),...
-                   0, 1, 10, 7, 0, 'Constituent'}, t, 0.2);
+                   0, 1, 10, 7, 0, 'Constituent', 0}, t, 0.2);
 T_mov.solve_tern_frap();
 norm_fac = 1/sum(diff(T_mov.x)...
                  .*(T_mov.sol(1, 1:end-1)+T_mov.sol(1, 2:end))/2);
@@ -400,31 +400,37 @@ end
 %% %%%%%%%%%%%%%%%%%%%%%% FRAP JUMP LENGTH DISTRIBUTION %%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 params = {-5, b(7/3, 10^-6), 0.5, e(7/3), 0, 1, 10, 7, 0, 'Constituent', 0};
-t = linspace(0, 5, 51);
+t = linspace(0, 2, 101);
 bp = 5;
 direc = -2; % 1: jump from left to right, -1: jump from right to left.
 
 %% 
-[ls, p_i_i] = frap_distro(t, params, 2, bp);
-[~, p_o_o] = frap_distro(t, params, -2, bp);
-[~, p_i_o] = frap_distro(t, params, 1, bp);
-[~, p_o_i] = frap_distro(t, params, -1, bp);
+[ls, p_i_i] = frap_distro(t, params, 2, bp, 1);
+[~, p_o_o] = frap_distro(t, params, -2, bp, 1);
+[~, p_i_o] = frap_distro(t, params, 1, bp, 1);
+[~, p_o_i] = frap_distro(t, params, -1, bp, 1);
+%%
+p_i_i(isnan(p_i_i)) = 0;
+p_o_o(isnan(p_o_o)) = 0;
+p_i_o(isnan(p_i_o)) = 0;
+p_o_i(isnan(p_o_i)) = 0;
 %%
 figure; hold on;
 plot(ls, sum(p_i_i, 2));
 plot(ls, sum(p_o_o, 2));
 plot(ls, sum(p_i_o, 2));
 plot(ls, sum(p_o_i, 2));
+save('/data/biophys/hubatsch/MatlabWorkspaces/prob_laplace_190121.mat');
 %% Save time points for comparison with python
 csvwrite('lb_in_out.csv', ls)
 csvwrite('prob_in_out.csv', p);
 %%
-KL = zeros(1, 20);
-for i = 1:20
+KL = zeros(1, 100);
+for i = 1:68
 i_t = 1+i;
 j_t = 2+i;
 integrand = @(p, q) abs(p.*log(p./q));
-integrate = @(x, y) sum(diff(x)'.*(y(2:end)+y(1:end-1))/2);
+integrate = @(x, y) nansum(diff(x)'.*(y(2:end)+y(1:end-1))/2);
 n_i = @(p1, p2, p3, p4) abs(integrate(ls, p1) + integrate(ls, p2)+...
                             integrate(ls, p3) + integrate(ls, p4));
 N_i = n_i(p_i_i(:, i_t), p_o_o(:, i_t), p_i_o(:, i_t), p_o_i(:, i_t));
@@ -435,18 +441,22 @@ KL(i) = integrate(ls, integrand(p_i_i(:, i_t)/N_i, p_i_i(:, j_t)/N_j))+...
      integrate(ls, integrand(p_i_o(:, i_t)/N_i, p_o_i(:, j_t)/N_j))+...
      integrate(ls, integrand(p_o_i(:, i_t)/N_i, p_i_o(:, j_t)/N_j));
 end
-plot(t(2:21), abs(KL));
-ylim([0, 1])
+plot(t(2:101), abs(KL)/0.02);
+% ylim([0, 1])
 %%
 figure; hold on;
-plot(ls, p_i_i(:, i_t)/N_i);
-plot(ls, p_o_o(:, i_t)/N_i);
-plot(ls, p_i_o(:, i_t)/N_i);
-plot(ls, p_o_i(:, i_t)/N_i);
+plot(ls, p_i_i(:, j_t)/N_j, 'b');
+plot(ls, p_o_o(:, j_t)/N_j, 'b');
+plot(ls, p_i_o(:, j_t)/N_j, 'b');
+plot(ls, p_o_i(:, j_t)/N_j, 'b');
+%%
+plot(T_mov.t, s_dot)
+hold on; xlim([0, 2]);
+plot(t(2:101), 0.1908/0.1885*0.61*abs(KL)/0.02);
 %% %%%%%%%%%%%%%%%%%%% INTEGRATION FUNCTION DEFINITIONS %%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-function [ls, p] = frap_distro(t, params, direc, bp)
+function [ls, p] = frap_distro(t, params, direc, bp, ind_delta)
 x0 = sort(5-abs(direc)/direc*(0:0.002:4.01));
 
 % Run simulations for 'delta' IC across outside
@@ -463,12 +473,12 @@ end
 T_mov = Ternary_model(0, 'FRAP', params, t, 0.2);
 T_mov.solve_tern_frap();
 ls = -abs(direc)/direc*(0.001:0.04:4);
-n_T=45;
+n_T=70;
 p = nan(length(ls), n_T);
 for j = 1:n_T
     tic
     parfor i = 1:length(ls)
-        p(i, j) = int_prob(ls(i), F, x0, direc, j, bp, 5, T_mov);
+        p(i, j) = int_prob(ls(i), F, x0, direc, j, bp, ind_delta, T_mov);
     end
     toc
 end