
A = zeros(50,50);
A(1,[1 2 3]) = [12 -2 1];
A(2,[1 2 3 4]) = [-2 12 -2 1];
A(49,[47 48 49 50]) = [1 -2 12 -2];
A(50,[48 49 50]) = [1 -2 12];
for i = 3:48
    A(i,[i-2 i-1 i i+1 i+2]) = [1 -2 12 -2 1];
end
b = ones(50,1)*5;

%% Solve using matrix inversion: 
AI = inv(A);
x1 = AI*b;

%% Solve using Mathews and Fink upper triangulation and back-substitution
%% function: 
x2 = uptrbk(A,b);

%% Solve using the Mathews and Fink Gauss-Seidel program

P = zeros(50,1);
delta = 1e-4*ones(50,1);
max1 = 10;

x3=gseid(A,b,P,delta, max1)


figure(1)
clf
spy(A)
title('sparsity pattern of A\n');

figure(2)
clf
stem(x1);
hold on
plot(x2, 'rd', 'MarkerSize', 12);
plot(x3, 'k');
title('solutions to Ax = b');


