A = magic(9);

b = [910; 1034; 1113; 1264; 1172; 981; 1060; 941; 750];

x_exact = A\b;

x = lu_solver(A,b);

dif = norm(x-x_exact);