n = 9

%% Set up Vandermode problem
x = (2:n+1)'
V = single(vander(x))

%% Random RHS
y = single(rand(n,1))

%% Compute the result using naive GE
%a = gauss_elim_vector(V, y)
a = gauss_elim_outer(V, y)
V*a-y

%% Plot the points and the curve
figure
plot(x,y,'o')
hold on
xs = linspace(1,n+2,201);
plot(xs,polyval(a, xs))
axis([1 n+2 -1 2])

%% Compute using backslash
a = V \ y
V*a-y