%computes expected times for random walk on 
%1,2,...,n+2 to reach 1 or n+2 for the first time

n = 5;

Q = zeros(n,n);      %transition matrix for random walk, cut off at edges
for i=1:n
    for j=1:n
        if abs(i-j) == 1
            Q(i,j) = -1/2;
        end
        if i == j
            Q(i,j) = 1;
        end
    end
end

A = [Q,ones(n,1)];   %augment A with all ones vector

%now we solve the tridiagonal system given by A

%forward elimination

for i=1:n-1
    r = A(i+1,i)/A(i,i)
    A(i+1,i) = 0;
    A(i+1,i+1) = A(i+1,i+1)-r*A(i,i+1)
    A(i+1,n+1) = A(i+1,n+1)-r*A(i,n+1)    %update far right column
end

%backward elimination

for i=n:-1:2
    r = A(i-1,i)/A(i,i);
    A(i-1,i) = 0;
    A(i-1,n+1) = A(i-1,n+1)-r*A(i,n+1);   %update far right column
end

%scale rows

for i=1:n
    r = A(i,i);
    A(i,i) = 1;
    A(i,n+1) = A(i,n+1)/r;
end

x = A(:,n+1);      %x(i) = expected time to reach 1 or n+2, starting at i+1
z = linspace(2,n+1,n);
plot(z,x)
ylabel('average first time to hit 1 or 7')
xlabel('starting point')