
clear all;
close all;

%c=.75;
c=1;
delx=.1;
a=-1;
b=5;
x=a:delx:b;
J=(b-a)/delx; % set number of spatial steps. The nr. of dividing points is J+1
delt=.1;
r=c*delt/delx; % c/rho.
nsteps=20;

u=f(x); % initialize the vector u.

% Forward difference
figure
for n=1:nsteps
    v=u;
    for j=1:J
        u(j)=(1+r)*v(j)-r*v(j+1);        
    end;
    if(n==1)
        plot(x,u,'g-+');
    end;
    hold on;
    if(n==nsteps)
        plot(x,u,'r-s')   
        hold on
        plot(x,sol(x,delt*nsteps)); % plot the solution with the method of characteristics  
        %xlim([0 2])
        %ylim([-1 1])
    end;
end;
title('The solution using a forward difference approx.')
xlabel('x','FontSize',14);
ylabel('u_{fwd}(x,t)','FontSize',14);
legend('n=1','n=nsteps','u_{exact}(x,4)')
hold off;

%break;

u=f(x);

% Upwind scheme
figure
for n=1:nsteps
    v=u;
    for j=2:J+1
        u(j)=(1-r)*v(j)+r*v(j-1);        
    end;
    if(n==1)
       plot(x,u,'g-+');
    end;
    hold on;
    if(n==nsteps)
        plot(x,u,'r-s') 
        hold on
        plot(x,sol(x,delt*nsteps)); % plot the solution with the method of characteristics
    end;
end;
title('The solution using an upwind difference approx.')
xlabel('x','FontSize',14);
ylabel('u_{uwd}(x,t)','FontSize',14);
legend('n=1','n=nsteps','u_{exact}(x,4)')
hold off;

figure
plot(x,f(x),'r'); % plot the initial solution
hold on
plot(x,sol(x,4)); % plot the solution with the method of characteristics
hold off;
title('Initial solution and Char.solution u(x(t),t)'); 
xlabel('x','FontSize',14);
ylabel('u(x,t)','FontSize',14);
legend('f','u_{char}(x,4)')

