% Script for solving Burger's equation
%  du/dt + u*du/dx = a*d2u/dx2, 0 < x < 1
% with periodic BCs using explicit RK4 time differencing and centered space
% differencing.

  L = 1;
  T=0.5;
  tplt = 0.1;  % Time between plots
  a = 0.01;
  nx = 100;
  dx = L/nx;
  x = (1:nx)*dx;

  u = 1+sin(2*pi*x);

  nplt = round(T/tplt);  
  cmax = max(u);  % Max wave propagation speed
  
  srmax = 2;  % Should be less than 2.82 RK4 stability limit for decay
  simax = 2;  % Should be less than 2.82 RK4 stability limit for oscillations

  % Need dt < min(dtr, dti)
  
  dtr = srmax*dx^2/(4*a);
  dti = simax*dx/cmax;
  dt = min(dti,dtr);

  % Adjust dt down so an integer number of timesteps fit into plot
  % time T

  ntplt = ceil(tplt/dt);
  dt = tplt/ntplt;
  
  subplot(2,2,1)
  plot(x,u,'r')
  xlabel('x')
  xlim([0 L])
  ylabel('u')
  hold on
  
  F = @dudt_burgers;
  
  % RK4 timestepping of du/dt = F
  
  for iplt = 1:nplt
  for n = 1:ntplt
    q1 = dt*F(u,x,a);
    q2 = dt*F(u+0.5*q1,x,a);
    q3 = dt*F(u+0.5*q2,x,a);
    q4 = dt*F(u+q3,x,a);
    u = u + (q1+2*q2+2*q3+q4)/6;
  end
  plot(x,u,'b')
  end

  hold off
