using SymPy

# Create an n x n (default 3), i,j elimination matrix
function E(i,j,n=3)
    E = sympy"eye"(n)           # start with the identity
    E[i,j] = - symbols("m$i$j")  # insert the negative multiplier
    E
end

# Create a symbolic 3x3 A matrix
  A = [symbols("A$i$j") for i=1:3, j=1:3]

E(2,1) 

using Interact

@manipulate for i=1:3, j=1:3
     E(i,j)
end

E(2,1)

inv(E(2,1))

E(2,1) * A   

E(2,1)^2 * A  

E(3,2) * A

E(2,1) * E(3,2)  * A

E(2,1) * A

E(3,2) * E(2,1) * A

E(3,2) * E(2,1)

E(2,1) * E(3,2)

E(2,1,5)

E1 = E(2,1,5) * E(3,1,5) * E(4,1,5)* E(5,1,5)

E(3,1,5) * E(2,1,5) * E(5,1,5)* E(4,1,5)  # Why does the order not matter

E1 # subtracts multiples of the first row

inv(E1) # INVERSE means add back the same multiples of the first row

E(3,2,5)

inv(E(3,2,5))

E1 * E(3,2,5)

inv(E1) * inv(E(3,2,5))

E1 * E(3,2,5)  # Why does this have a simple looking answer?

inv(E1) * inv(E(3,2,5)) # Why does this have an even slightly simpler looking answer?

 E(3,2,5) * E1 # Why doesn't this have a simple looking answer?

E2 = prod( E(i,2,5) for i=3:5)

E3 = prod( E(i,3,5) for i=4:5)

E4 = E(5,4,5)

E1 * E2 * E3 * E4 # Why is this simple?

L = inv(E1) * inv(E2)  * inv(E3) * inv(E4) # Why is this simple?  This is the L Matrix!!

E4 * E3 * E2 * E1 # Why is this a mess?  Good thing we don't need it in Gaussian Elimination

inv(L)  # Right this is the inv(L) , not how the inverse of a triangular matrix isn't so pretty in general

L

b(i) = symbols("b$i")

[b(i) for i=1:5]

y = L\[b(i) for i=1:5]

y[1] = b(1)

y[2] = b(2) - L[2,1]*y[1]

y[3] = b(3) - L[3,1]*y[1] - L[3,2]*y[2]

y = [b(i) for i=1:5]
for i = 1:5, j=1:i-1  
    y[i] -= L[i,j]*y[j]
end

y