title: How to do Gaussian elimination in R (do not use "solve")
body: <pre><code>A &lt;- matrix(c(2,-5,4,1,-2.5,1,1,-4,6),byrow=T,nrow=3,ncol=3)
  b &lt;- matrix(c(-3,5,10),nrow=3,ncol=1)
  p &lt;- nrow(A)
  U.pls &lt;- cbind(A,b)
  for (i in 1:p){
      for (j in (i+1):(p+1)) U.pls[i,j] &lt;- U.pls[i,j]/U.pls[i,i]
      U.pls[i,i] &lt;- 1
      if (i &lt; p) {
         for (k in (i+1):p) U.pls[k,] &lt;- 
                     U.pls[k,] - U.pls[k,i]/U.pls[i,i]*U.pls[i,]
      }
  }
  U.pls
  x &lt;- rep(0,p)
  for (i in p:1){
    if (i &lt; p){
       temp &lt;- 0
       for (j in (i+1):p) temp &lt;- temp + U.pls[i,j]*x[j]
       x[i] &lt;- U.pls[i,p+1] - temp
    }
    else x[i] &lt;- U.pls[i,p+1]
  }
  x
  </code></pre>
  
  <h1>output</h1>
  
  <pre><code>&gt; U.pls
       [,1] [,2] [,3] [,4]
  [1,]    1 -2.5    2 -1.5
  [2,]    0  1.0 -Inf  Inf
  [3,]    0  0.0    1  NaN
  &gt; x
  [1] NaN NaN NaN
  </code></pre>
  
  <p>In this way, I cannot solve solution in some case. Even though I know the reason why the error happens mathematically, I cannot fix the error in R.
  Help me with some fix. Thank you in advance.</p>