save_figures = 0;

disp('Graphical interpretation:');
dx = L_x:0.01:R_x;
figure;

% plotting initial condition
plot(dx, IC(dx), 'LineWidth', 1, 'Color', 'blue', 'LineStyle', '--'); hold on;

% plotting exact solution
exact_solution = IC(dx - a_advec * T); % exact solution to for plotting
plot(dx, exact_solution, 'LineWidth', 1, 'Color', 'black', 'LineStyle', '--'); hold on;

% the correct way of plotting computed solution using nodal representation
dr_plot = -1:0.01:1; 
LagPols_cell = LagrangePolynomialBasis(r);
LagPols_values_plot = zeros(N_p,length(dr_plot));
for i = 1:N_p
    LagPols_values_plot(i,:) = polyval(LagPols_cell{i},dr_plot);
end
for i = 1:K
    element = E(:,i); x_l_ID = element(1); x_r_ID = element(2);
    x_l = P(x_l_ID); x_r = P(x_r_ID);
    x_plot = x_l + ((1 + dr_plot) / 2) * (x_r - x_l);
    coef = u(:,i)';
    LagPol_values_plot = coef * LagPols_values_plot;
    plot(x_plot, LagPol_values_plot, 'LineWidth', 1.5, 'Color', 'red'); hold on;
end

% figure annotation
axis(axis_limits); xlabel('x'); ylabel('u');
legend('IC','exact','computed','Location','NorthWest');

%if alpha_flux == 1 % pro central flux
%    title_string = ['solution in time t = T = ',num2str(T),'s: computed = red, exact = black; "central flux" chosen'];
%end
%if alpha_flux == 0 % pro upwind flux
%    title_string = ['solution in time t = T = ',num2str(T),'s: computed = red, exact = black; "upwind flux" chosen'];
%end
%title(title_string);

grid on;

% saving figure into a file
if save_figures == 1
    if alpha_flux == 1 % for central flux
        filename_string = ['Advection1D_solution_for_K=',num2str(K),',_N=',num2str(N),'_and_central_flux'];
    end
    if alpha_flux == 0 % for upwind flux
        filename_string = ['Advection1D_solution_for_K=',num2str(K),',_N=',num2str(N),'_and_upwind_flux'];
    end
    saveas(gcf,filename_string,'pdf');
    saveas(gcf,filename_string,'png');
end