% Script for solving Burger's equation 
%  du/dt + u*du/dx = -a*d2u/dx2, 0 < x < 1
% with periodic BCs using splitting, with RK4 time differencing/centered space
% differencing for the advection term and Crank-Nicolson for the
% diffusion term.

  L = 1;
  T=0.5;
  tplt = 0.1;  % Time between plots

  a = 0.01;
  nx = 100;
  dx = L/nx;
  x = (1:nx)'*dx;    % Make x and u column vectors for convenience with C-N
  u = 1+sin(2*pi*x);

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

  dti = simax*dx/cmax;

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

  ntplt = ceil(tplt/dti);
  dt = tplt/ntplt;
  
  subplot(2,2,1)
  plot(x,u,'r')
  xlabel('x')
  xlim([0 L])
  ylabel('u')
  hold on
  
  for iplt = 1:nplt
    for n = 1:ntplt
  
      % RK4 timestepping of du/dt = -u*du/dx
  
      q1 = dt*dudt_burgers_adv(u,x);
      q2 = dt*dudt_burgers_adv(u+0.5*q1,x);
      q3 = dt*dudt_burgers_adv(u+0.5*q2,x);
      q4 = dt*dudt_burgers_adv(u+q3,x);
      u = u + (q1+2*q2+2*q3+q4)/6;
    
      % Crank-Nicolson on du/dt = a*du2/dx2

      nu = a*dt/dx^2;
      jp1 = [2:nx 1];
      jm1 = [nx 1:(nx-1)];
      e = ones(nx,1);
      M = spdiags([-0.5*nu*e e*(1+nu) -0.5*nu*e],-1:1,nx,nx);
      M(1,nx) = -0.5*nu;
      M(nx,1) = -0.5*nu;
    
      rhs = u*(1-nu) + 0.5*nu*(u(jp1)+u(jm1));
      u = M\rhs;
              
  end
  plot(x,u,'b')
  end

  hold off
