#set of inputs
E1 <- function (x0, x1, x2) 0.3 * x0 + -0.2 * x1 + 10 * x2 + -71.4;
E2 <- function (x0, x1, x2) 3 * x0 + -0.2 * x2 + -0.1 * x1 + -7.85;
E3 <- function (x0, x1, x2) 0.1 * x0 + 7 * x1 + -0.3 * x2 + 19.3;
system <- list(E1, E2, E3);

CreateCoeffMatrix <- function (input) { #create a blank coefficient matrix with rownames
  coeffMatrix = matrix (0, nrow=length(input), ncol=length(input)+1, byrow=TRUE)
  rownames = c()
  for(i in 1:nrow(coeffMatrix)) {
    rownames[i] = i
  }
  rownames(coeffMatrix) <- rownames
  return(coeffMatrix)
}

AugCoeffMatrix <- function (system) { #create the augmented coefficient matrix
  augcoeffmatrix = CreateCoeffMatrix(system)
  variables = c() #initially empty list of variables

  for (i in 1:length(system)) { #parse the variables
    ex = deparse(system[i])[1]
    ex = gsub("list", "", ex) #remove excess characters
    ex = gsub("function", "", ex) #remove excess characters
    ex = gsub("[\\(\\) ]", "", ex, ",") #remove excess characters
    ex = strsplit(ex, ",") #split variables by commas
    varlist = unlist(ex) #store in varlist
    for (j in 1:length(varlist)) { #check variables if already in variable list. if not, store in variables
      if (varlist[j] %in% variables == FALSE) {
        variables = c(variables, varlist[j])
      }
    }
  }

  if (length(variables) != length(system)){ #return NA if the number of variables is not equal to the number of equations
    return (NA)
  }

  colnames = variables
  colnames = c(colnames, "RHS")
  colnames(augcoeffmatrix) <- colnames #provide column names for the augmented coefficient matrix

  vals <- list()
  eqList = c()

  for (i in 1:length(system)) { #parse expressions
    eq = deparse(system[i])[2]
    eqList = c(eqList, eq) #copy of equations
    eq = gsub("\\*", "", eq) #remove excess characters
    eq = gsub("\\s", "", eq) #remove excess characters
    eq = gsub("\\(", "", eq) #remove excess characters
    eq = gsub("\\)", "", eq) #remove excess characters
    eq = gsub("x[0-9]", "", eq) #remove excess characters
    eq = strsplit(eq, "\\+") #split by the + sign
    vals = c(vals, eq) #store numbers in vals
  }
  for (i in 1:nrow(augcoeffmatrix)) { #store coefficients in augcoeffmatrix
    nums = as.numeric(unlist(vals[[i]])) #gets each coefficient then turn it into a numeric typing
    numReference = unlist(vals[[i]]) #gets each coefficient as a string for string manipulation
    for (j in 1:length(variables)) {
      locVar = regexpr(numReference[j], eqList[i]) #gets location of coefficient in the expression
      varLen = attr(locVar, "match.length") #gets the length of the coefficient
      varEx = substr(eqList[i], locVar[1]+varLen+3, locVar[1]+varLen+4) #gets the corresponding variable to the coefficient

      for(k in 1:ncol(augcoeffmatrix)){ #store in augcoeffmatrix
        if (varEx == colnames[k]) { #if the coefficient's variable is the same as the matrix column, store
          augcoeffmatrix[i,k] = nums[j]
        }
        if (k == ncol(augcoeffmatrix)) { #if the value is the last element, it is the constant. negate and add to the RHS column
          augcoeffmatrix[i,k] = -nums[k]
        }
      }
    }
  }
 returnList = list(augcoeffmatrix, variables) #lists the return values
 return(returnList)
}