In this paper, a new method is proposed to find the approximate solutions to fully fuzzy systems of linear equations (FFSLEs). The technique integrates a bisection method with Fuzzy Numerical Simulation (FNS). The procedure starts with generating single values of fuzzy parameters and solving the resulting crisp problems repeatedly to determine the lower and upper bounds of the solutions. After computing the mean lower and upper bound values, the obtained supremum and infimum values are considered to be the lower and upper bounds of the solutions, respectively. It is attempted to improve solutions by considering an error function related to the sum of the absolute differences between the corresponding lower and upper bounds of the left and right sides of the equalities. When very large intervals are obtained for the solutions, the bisection algorithm is applied to reduce the error value. The method intends to solve square systems of large dimensions for arbitrary fuzzy numbers (FNs) by removing non-negativity confinements of the variables and/or coefficients to be more realistic. After the computational method is presented thoroughly, some benchmark examples are finally provided.