#Nathan Fox
#Homework 26
#I give permission for this work to be posted online

#Read procedures from class
read(`C26.txt`):

Help:=proc():
    print(` ExpD2(M,N,pleft:=1/4,pright:=1/4,pdown:=1/4,pup:=1/4) `):
end:

##PROBLEM 2##

#ExpD2(M,N,pleft:=1/4,pright:=1/4,pdown:=1/4,pup:=1/4):
#inputs positive integers M and N, and outputs a list
#of lists, let's call it L, such that
#L[i][j] (0 ≤ i ≤ M-1   ,   0 ≤ j ≤ N-1 )
#is the expected number of unit steps left for a random walker
#who starts at [i,j] and moves, in a unit step, with probabilities [pleft, pright, pdown, pup] to one of its immediate neighbors
#[i-1,j],[i+1,j],[i,j-1],[i,j+1] (repsectively)
#and ends the walk when reaching one of the edges
#i=0, i=M, j=0, j=N
#(Note: if pleft+pright+pdown+pup<>1, then they are all scaled
#appropriately.  Also, they all default to 1/4
ExpD2:=proc(M, N, pleft:=1/4, pright:=1/4, pdown:=1/4, pup:=1/4)
local P, eq, var, var1, f, i, j:
    P:=[pleft, pright, pdown, pup]/(pleft+pright+pdown+pup):
    var:={seq(seq(f[i,j], i=1..M-1), j=1..N-1)}:
    eq:={seq(seq(f[i,j]=pleft*f[i-1,j] + pright*f[i+1,j] +
        pdown*f[i,j-1] + pup*f[i,j+1] + 1, i=1..M-1), j=1..N-1)}:
    eq:=subs({seq(f[0,j]=0, j=1..N-1), seq(f[M,j]=0, j=1..N-1),
        seq(f[i,0]=0, i=1..M-1), seq(f[i,N]=0, i=1..M-1)}, eq):
    var:=solve(eq, var):
    return [seq([seq(subs(var, f[i,j]), j=1..N-1)], i=1..M-1)]:
end: