function method_of_lines

%% Output time
final_time = 5;
%% Deltat
deltat = 0.01;
CFL = 0.75;
%% Number of nodes in space
n = 5000;
%%Compute deltax
x_min = 0.0;
x_max = 1.0;
deltax = (x_max -x_min)/n;
[x,rho] = set_initial_condition(x_min ,deltax ,n);
current_time = 0.0;
drhodt = zeros(n,1);
y = zeros(n,1)

    while current_time < final_time
        deltat = find_deltat(rho , deltax , CFL);
        if (current_time+deltat > final_time)
            deltat = final_time - current_time;
        end

        current_time = current_time + deltat;
        y(1) = 1 - 1*exp(-8*exp(-3*rho(1)));

        for i = 2:n-1
            y(i) = 1 - 1*exp(-8*exp(-3*rho(i)));
        end

        y(n) = 1 - 1*exp(-8*exp(-3*rho(n)));

        drhodt(1) = -1/(deltax) * (y(2) - y(1));
        for i = 2:n-1
            drhodt(i) = -1/(deltax) * (y(i+1) - y(i));
        end
        drhodt(n) = -1/(deltax) * (y(1) - y(n));

        rho = rho + deltat*drhodt;

        figure (1);
        plot(x,rho ,'-o', 'LineWidth', 5);
        time = strcat(num2str(current_time), 's');
        tex = strcat('t = ', time);
        text(0.6,2, tex);
        legend('Finite -Difference')
        title('Time marching a P.D.E.')
        xlabel('x')
        ylim ([0 ,2.5]);
        ylabel('Traffic Density')
        grid on
        drawnow
    end
        figure (2);
        plot(x,y,'-o', 'LineWidth', 5);
        legend('Finite -Difference')
        title('Time marching a P.D.E.')
        xlabel('x')
        ylim ([0 ,1.5]);
        ylabel('Traffic Velocity')
        grid on
        drawnow

end
function [x,rho] = set_initial_condition(x_min ,deltax ,n)
x = zeros(n,1);
rho = zeros(n,1);
for i = 1:n
    x(i) = x_min + (i)*deltax;
    rho(i) = 1;
end
for i = 1:n
    if 0.25 <= x(i) && x(i) <= 0.75
        rho(i) = -16*x(i)^2 + 16*x(i) - 2.0;
    end
end
end
function deltat = find_deltat(rho ,deltax ,CFL)
n = size(rho ,1);
deltat = 1E100;
%find minimum dt
for i = 1:n
    max_speed = abs(-24*exp(-3*rho(i) -8*exp(-3*rho(i))));
    deltat = min(CFL*deltax/max_speed , deltat);
end
end