Implementing a markov chain for a probability dependent random walk

Consider a particle hopping on a one-dimensional lattice with lattice spacing s. The particle starts at the lattice site x0. The particle jumps left or right depending on a position dependent probability distribution. Concretely, I have been able to determine the following recursion relation for the particle motion

For a given N we have .

Here is the probability that the particle reaches position after a total of N jumps. Φ is the normal cumulative distribution function, with mean μ and standard deviation σ. is the Kronecker delta symbol.

To avoid confusion: in my implementation I renamed position to velocity.

My implementation seems to be working fine, but the one-indexing and array-indexing is throwing me off.

If I run an animation with the particle moves to the right towards the mean, and vice versa for . However, if I look at the output matrix called "distribution" I notice that all the elements are zero above column 158 and below column 55. Of course this might be due to rounding in matlab, but I became uncertain if I somehow messed up the indices in my recursion relation.

Unfortunately, my question then reduces to the mundane: are my implemented indices correct in the following loop?

for N = 2:Nmax
    for n = n0-(N-1) :2: n0+(N-1)
        distribution(N,n)=(1-kronDel(n0-(N-1),n))*distribution(N-1,n-1)*(1-normcdf((V(n-1)-mean)/sigma))...
                          +(1-kronDel(n0+(N-1),n))*distribution(N-1,n+1)*normcdf((V(n+1)-mean)/sigma);
    end

end

Entire code:

%% Settings
clc
close all
clear all
set(0,'defaulttextinterpreter','latex')
set(0,'defaultAxesTickLabelInterpreter','latex');
set(0,'defaultLegendInterpreter','latex');
set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultAxesFontSize',15)

N=100;
mean = 400;
sigma = 30;
s = 5;
V0 = 370;

[distribution, Velocity]=f(N,mean,sigma,s,V0);

figure(1)
for i = 1:N
    plot(Velocity,distribution(i,:))
    title(i)
    xlabel('Velocity [m/s]')
    ylabel('Probability')
    pause(0.1)
    drawnow
end

function [distribution_up_down_markov, Velocity] = f(Nmax,mean,sigma,s,V0)

%N=1 (0) case.
%Possible velocities
V = V0-s*Nmax:s:V0+s*Nmax;
%V(Nmax+1)=V0;
distribution = zeros(Nmax,length(V));
distribution(1,Nmax+1)=1;
n0 = Nmax + 1;

for N = 2:Nmax
    for n = n0-(N-1) :2: n0+(N-1)
        distribution(N,n)=(1-kronDel(n0-(N-1),n))*distribution(N-1,n-1)*(1-normcdf((V(n-1)-mean)/sigma))...
                          +(1-kronDel(n0+(N-1),n))*distribution(N-1,n+1)*normcdf((V(n+1)-mean)/sigma);
    end

end


distribution_up_down_markov=distribution;
Velocity = V;

end

function d=kronDel(j,k)

d=j==k;

end

Answer 1

[VBBV]

for n = 1-n0 :2: n0-1

[VBBV]

From the given relation of N the inner loop index may be changed to take account of every 2 element in the range