Why is my Gauss-Seidel function slightly off?

I've written a Gauss Seidel function that takes in a coefficient matrix, a vector array holding the constants, and an integer n, where n-1 is the length of the matrix. It returns an array that holds the variable solutions. I've found that my variables are slightly off. This is one of my test cases:

double matrix[5][5] = 
{ {-1.601675, 0.700000, 0.000000, 0.000000, 0.000000},
  {0.700000, -1.201675, 0.500000, 0.000000, 0.000000},
  {0.000000, 0.500000, -0.801675, 0.300000, 0.000000},
  {0.000000, 0.000000, 0.300000, -0.401675, 0.100000},
  {0.000000, 0.000000, 0.000000, 1.000000, -1.008375},
};

double sol[5] = 
{
  -180.041874,
  -0.041874,
  -0.041874,
  -0.041874,
  -0.209372,
};

I expect to get (and ran this system of equations through Wolfram Alpha) the following for my unknowns:

double T_expected[5] = 
{
  198.298282,
  196.616925,
  194.965693,
  193.373744,
  192.046373,
};

However, running the equations through my function, I get these values instead.

double T[5] = 
{
  199.674915,
  199.257945,
  198.676137,
  197.711841,
  196.277411,
};

Here is my function:

#define EPSILON 1e-06
double *runGaussSeidel(double **matrix, double *sol, int n)
{
    int i, j, converged, m = n - 1;
    double sum1 = 0.0, sum2 = 0.0, temp = 0.0;
    double *T = malloc(sizeof(double) * m);

    // Initial guess.
    for (i = 0; i < m; i++)
        T[i] = 0.0;

    do {   
        converged = 1; 
        for (i = 0; i < m; i++)
        {
            sum1 = sum2 = 0.0;
            for (j = 0; j < i; j++)
            {
                sum1 += matrix[i][j] * T[j];
            }   
            for (j = i + 1; j < m; j++)
            {
                sum2 += matrix[i][j] * T[i];
            }

            double temp = T[i];

            T[i] = (sol[i] - sum1 - sum2) / matrix[i][i];                                                                   
        

            if (T[i] - temp >= EPSILON)
            {
                converged = 0;
            }                    
        }                                                                      
    } while (!converged); 
    
    return T;
}

I think you have a bug - you change T[i] and then use it for the next variable/iteration. You have to save previous solution vector before any assignments. Also why n-1 I noted earlier.

Nevertheless, your T_expected[5] doesn't look right as well. Did you try to test it?

Here is my code,

#include <alloca.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

void iteration(
    const int n, 
    double (*mtx)[n][n],
    const double* rhs,
    double* res
) {
    // save solution values to use in iteration
    double* tr = alloca(n * sizeof(double));
    memcpy(tr, res, n*sizeof(double));
    
    // actual iteration
    for(int k = 0; k != n; ++k) {
        double s = 0.0;
        for(int i = 0; i != n; ++i) {
            s += (*mtx)[k][i] * tr[i];
        }
        double diff = (rhs[k] - s)/(*mtx)[k][k];
        res[k]     += diff;
    }
}

double test(
    const int n,
    double (*mtx)[n][n],
    const double* rhs,
    const double* res
) {
    double d = 0.0;
    
    for(int k = 0; k != n; ++k) {
        double s = 0.0;
        for(int i = 0; i != n; ++i) {
            s += (*mtx)[k][i]*res[i];
        }
        
        d += fabs(rhs[k] - s);
    }
    return d;
}

int main(void) {
    double matrix[5][5] = {
        {-1.601675, 0.700000, 0.000000, 0.000000, 0.000000},
        {0.700000, -1.201675, 0.500000, 0.000000, 0.000000},
        {0.000000, 0.500000, -0.801675, 0.300000, 0.000000},
        {0.000000, 0.000000, 0.300000, -0.401675, 0.100000},
        {0.000000, 0.000000, 0.000000, 1.000000, -1.008375}
    };

    double sol[5] = {
        -180.041874,
        -0.041874,
        -0.041874,
        -0.041874,
        -0.209372
    };

    double res[5] = { 0, 0, 0, 0, 0};
    
    for(int k = 0; k != 1000; ++k) {
        iteration(5, &matrix, sol, res);
        printf("%f\t%f\t%f\t%f\t%f\n", res[0], res[1], res[2], res[3], res[4]);
    }

    printf("%g\n", test(5, &matrix, sol, res));

    return 0;
}

and after 1000 iterations solution is

198.567523  197.141091  195.720761  194.306738  192.900568

it computes sum of absolute differences (test() function) being equal to

7.65291e-14