% Here we solve u_t+au_x=0, u(x,0)=sin(2*pi.*x)
% which has the exact solution: u(x,t)=sin(2*pi.*(x-at))
% Given the characteristics and the solution, we assume periodicity of U

clear;clf;
% Fix the time, number of space steps and mesh size; 
T=2; J=40; h=1/J; 
% Fix the mesh spacing and convection term
a=0.5;
% the CFL condition is satisfied
lambda=1; k=lambda*h;
% Make sure that an integer number of time steps are required.
N=T/k; N=round(N); k=T/N;
x=[0:h:1]';
t=[0:k:T]';
U=zeros([J+1,N+1]);
U(:,1)=sin(2*pi.*x);

tic
for n=1:N

% Upwind scheme The exact boundary condition is U(1,n+1)=sin(-2*pi*a*(n)*k);
% Note the sever loss in dissipation when compared with the Lax-Wendroff scheme.
U([2:J+1],n+1)=(1-a*lambda)*U([2:J+1],n)+a*lambda*U([1:J],n);
U(1,n+1)=U(J+1,n);

%% Lax-Wendroff
%U([2:J],n+1)=U([2:J],n)-a*lambda/2*(U([3:J+1],n)-U([1:J-1],n))+a^2.*lambda.^2/2*(U([3:J+1],n)-2*U([2:J],n)+U([1:J-1],n));
%U(1,n+1)=U(J+1,n);
%U(J+1,n+1)=U(J,n);

end
toc
% The exact answer

for n=1:N+1
V(:,n)=sin(2*pi.*(x(:)-a*(n-1)*k));
end

for j=1:N+1
 	   plot(x,U(:,j));axis([0,1,-1,1]);
       M(j) = getframe;
 	end
% 	movie(M)
 
%
%hold on; for i=1:81;plot(x,U(:,i)-V(:,i));end
% The numerical solution
%subplot(2,1,1)
%surf(x,t,U')
% The error
%subplot(2,1,2)
%surf(x,t,U'-V')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% A convection problem
%
% Used in 4H Numerical Analysis
%
% J.F.Blowey
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
