function p2
  A = [3.02,-1.05,2.53;4.33,0.56,-1.78;-0.83,-0.54,1.47];
  b = [-1.61;7.23;-3.38];

  disp('------------------------------------------------');
  disp('Part a) -- original matrix')
  disp('------------------------------------------------');
  sol1 = gauss_back(A, b);
  %sol_no_rounding = gauss_back_noround(A, b)
  exact_solution = A \ b

  disp('------------------------------------------------');
  disp('Part b) -- changed A(1,1)')
  disp('------------------------------------------------');
  A2 = A;
  A2(1,1) = 3.01;
  sol2 = gauss_back(A2, b);
  percentage_change = abs(100*(sol2-sol1)./sol1)
  %sol_no_rounding = gauss_back_noround(A2, b)
  exact_solution = A2 \ b

  disp('------------------------------------------------');
  disp('Part c) -- changed b(3)')
  disp('------------------------------------------------');
  b2 = b;
  b2(3) = -3.39
  sol2 = gauss_back(A, b2);
  percentage_change = abs(100*(sol2-sol1)./sol1)
  %sol_no_rounding = gauss_back_noround(A, b2)
  exact_solution = A \ b2
end
% -----------------------------------------------------------
function sol = gauss_back(A, b)
  Ab = [A,b]

  [n, dummy] = size(A);

  for step = 1:n-1
    disp('------------------------------------------------');
    disp(sprintf('Begin step %d', step));
    pivot = A(step, step)
    for i = step+1:n
      disp(sprintf('row(%d) <- row(%d) - %g/%g * row(%d)', i, i, A(i,step), pivot, step));
      factor = myround(A(i,step)/pivot);
      A(i,step:n) = A(i,step:n) - factor*A(step,step:n);
      b(i) = b(i) - factor*b(step);
    end

    A = myround(A);
    b = myround(b);

    Ab = [A,b]
  end

  disp('------------------------------------------------');
  disp('Backsubstitution');
  sol = b;
  for step = 0:n-1
    row = n-step;
    sol(row) = sol(row) - A(row,(row+1):n)*sol(row+1:n);
    sol(row) = sol(row) / A(row,row);
    sol = myround(sol);
  end

  sol
  residual = A*sol-b
end
% -----------------------------------------------------------
function sol = gauss_back_noround(A, b)
  [n, dummy] = size(A);

  for step = 1:n-1
    pivot = A(step, step);
    for i = step+1:n
      factor = A(i,step)/pivot;
      A(i,step:n) = A(i,step:n) - factor*A(step,step:n);
      b(i) = b(i) - factor*b(step);
    end
  end

  sol = b;
  for step = 0:n-1
    row = n-step;
    sol(row) = sol(row) - A(row,(row+1):n)*sol(row+1:n);
    sol(row) = sol(row) / A(row,row);
  end
end
% -----------------------------------------------------------
function y = myround(x)
  y = round(x*100)/100;
end
