dt = 0.1;
dx = 1;

x = 0:dx:600;
nx = length(x);
nt = 1500;

c = 1;
D=1;

u = zeros(nx,1);
u(1:20) = 1;
u = u+1;

for t = 1:nt
    u2 = u;
    for i = 2:nx-1
        u(i) = u2(i) - c*dt/dx*(u2(i)-u2(i-1)) + D*dt/dx^2*(u2(i+1)-2*u2(i)+u2(i-1)); %c=U(i) ??
        %dU/dt=-c*dU/dx+2*dU^2/dx^2 %lin convection+diff
    end
    u(1) = 1;
    u(end) = 1;
    plot(u);
    ylim([0 3]);
    pause(0.01);
end