function AdvectionStable(a,dt,dx,max_time)

% Creates funtion input x = (0 dx 2*dx ... 10)^T
% dx = .05;
x = 0:dx:10;
Nu = a*dt/dx;
n_dim = length(x);
time_array = 0:dt:max_time;
time_index = length(time_array);

% Find IC

[IC] = AdvectionStableIC(dx);

plot(x,IC)
pause(.1)

Ones = linspace(1,1,n_dim-1);
A = eye(n_dim) - Nu^2*eye(n_dim) + (.5)*(Nu^2 + Nu)*diag(Ones,-1) + (.5)*(Nu^2 - Nu)*diag(Ones,-1);
A(1,n_dim)=.5*(Nu^2+Nu); 
A(n_dim,1)=.5*(Nu^2-Nu);

% Solved and plotted with Euler 

U = IC';

for i = 0:1:time_index;
    U = A*U;
    plot(x,U(1:n_dim));
    axis([0 10 -.1 2]);
    pause(.001)
    M(i+1) = getframe(gcf);
end    
% Remove Comment if you want to produce a movie.
% movie2avi(M,'AdvectionMovie.avi');


% For debugging
% dx
% dx
% a
% Nu
% A

end
