%
%  finite difference code:
%
%    u  + c(x) u  = 0
%     t         x   
%

a = 0.5;

% grid refinement, larger "work" means finer grid
work = 8;

%  grids (0 includes zero pt)

M = 16*work;    h = 2*pi/M;
N = 48*work;    k = 2*pi/N;

x  = h:h:2*pi;   x0 = [0 x];
t  = k:k:2*pi;   t0 = [0 t];

% initial value & c(x) vector

u  = sin(x + a*(1-cos(x)));  u0 = [0 u];
c  = 1./(1 + a*sin(x));      c0 = [1 c];

fc = (k/h)*c;

% solution matrix (begin with IVs)
u0_old = 0;
u_old  = u;

u_mn   = [0 u];

% iterate on t_n using finite difference formula
for n1 = 1:N
  u0_new  = sin(-n1*k);
  u_new   = u_old - fc.*(diff([u0_old u_old]));

  u_mn    = [u_mn; [u0_new u_new]];

  u0_old  = u0_new;
  u_old   = u_new;
end

% plot level curves of solution:  u_mn = u(x,t)

figure(1);               clf
subplot(2,1,1)
plot(x0,u_mn(1,:),'r')
hold on
plot(x0,u_mn(21:20:161,:),'k')
axis([0 2*pi -1.5 1.5])
title('\bf u(x,t) at time intervals of (2\rm\pi\bf/16)')
xlabel('\bf x');             ylabel('\bf u(x,t), t=0 is red')

subplot(2,1,2)
plot(x,fc)
axis([0 2*pi 0 2])           
title('\bf CFL number (k*c(x)/h)')
xlabel('\bf x');             ylabel('\bf CFL #  \rm(>1, unstable)')
axis([0 2*pi 0 2]);
hold on
plot([0 9],[1 1],'--')