% this m-file solves problem 6 from HW#2

clear all;
n = 7                                 % we solve the largest system here.
eps = 10^(-3)                         % with the given epsilon

h = 2^(-n)
N = 2^n -1

% the matrix A:
% we store the diagonal elements in the vector a, and the off diagonals in
% vectors b,c 
a = ones(1,N)*(2*eps/h^2+1)
b = ones(1,N)*(-eps/h^2)
c = ones(1,N)*(-eps/h^2)

% the right-hand side f:
x = h:h:1-h                            % create the coordinates
f = 2*x+1

% we now solve Au = f:

% first, forward sweep, which removes the lower diagonal
% and modifies the right hand side accordingly
for i = 2:N
    coeff = -b(i) / a(i-1);
    b(i) = 0;
    a(i) = a(i) + coeff * c(i-1);
    f(i) = f(i) + coeff * f(i-1);
end

% second, a backward sweep, which produces our solution u
u = ones(1,N);
u(N) = f(N) / a(N);
for i = N-1:-1:1
    u(i) = (f(i) - c(i)*u(i+1)) / a(i);
end

plot(x,u,'d','LineWidth',1.5,'Color',[1/n,1/n,2/n]);
axis([0,1,0,3]);