n = 5;
N = 2*n+1;    %random walk on {-n,-n+1,...,0,...,n-1,n}
q = 4/5;      %q = probability to move to the left on {-n,-n+1,...,-1}
              %q = probability to move to the right on {1,2,...,n}

P = zeros(N,N);   %transition probability matrix
for i=2:n
    P(i,i-1) = q;
    P(i,i+1) = 1-q;
end
P(n+1,n) = 1/2;
P(n+1,n+2) = 1/2;
for i=n+2:N-1
    P(i,i+1) = q;
    P(i,i-1) = 1-q;
end
P(1,2) = 1;
P(N,N-1) = 1;

mat = zeros(N,N-1);
x = linspace(-n,n,N);
y = linspace(-n,n,N);

for i=1:N
    Q = P;
    Q(i,:) = [];
    Q(:,i) = [];
    R = transpose(inv(eye(N-1)-Q)*ones(N-1,1));
    mat(:,i) = [R(1:i-1) 0 R(i:N-1)];
end

subplot(2,1,1)
z1 = linspace(-n,n-1,N-1);
z2 = linspace(-n+1,n,N-1);
u1 = zeros(1,N-1);
u2 = zeros(1,N-1);
for i=1:N-1
    u1(i) = P(i,i+1);
    u2(i) = P(i+1,i);
end
plot(z1,u1,'-b',z2,u2,'-r')
legend('probability to move right', 'probability to move left')
title('random walk transition probabilities')

subplot(2,1,2)
imagesc(x,y,mat)
set(gca,'YDir','normal')
colorbar
xlabel('x')
ylabel('y')
title('average time to reach y, starting at x')

    
    