function [resid, newu] = redblack(u, nx, h, dt, nguard, rhs, weight)
  
  odd = nguard+1:2:nx-nguard+1
  even = nguard:2:nx-nguard+1
 
  newu  = zeros(u); 
  du = zeros(u);
  resid = zeros(u);
  
  // red pass
  du(even) = 1./(h*h) * (u(even+1) -2.*u(even) + u(even-1));
  resid(even) = +rhs(even) + du(even);
  du(even) = dt*(du(even) + rhs(even));
  
  // black pass
  du(odd) = 1./(h*h) * (u(odd+1) -2.*u(odd) + u(odd-1));
  resid(odd) = +rhs(odd) + du(odd);
  du(odd) = dt*(du(odd) + rhs(odd));

  newu = u + weight*du; 

endfunction

function [source] = sourceterm(x)
  source = sin(4*x);
endfunction

// enforce  periodic boundary conditions
function [withbc] = setbc(withoutbc, nbc)
  n = max(size(withoutbc));
  withbc = withoutbc;
  for i=1:nbc 
      withbc(i) = withoutbc(n-2*nguard+i);
      withbc(n-i+1) = withoutbc(2*nguard-i+1);
  end
endfunction
    
// add bcs to something only containing the interior
function [withbc] = addbc(withoutbc,nbc)
  withbc = setbc([zeros(1,nbc), withoutbc, zeros(1,nbc)],nbc);
endfunction

// linear interpolation
function [ucoarse] = coarsen(u,nbc)
  [dmy,nu] = size(u);
  nc = (nu-1)/2;
  uc = zeros(dmy,nc)
  for i=2:nc-1
	uc(i) = (u(2*i-1)+2.*u(2*i)+u(2*i+1))/4.;
  end
  // strip old bc, add new bc
  uc = uc(nbc/2+1:nc-nbc/2);
  ucoarse = addbc(uc,nbc);
endfunction

//linear interpolation
function [ufine] = finer(u,nbc)
  [dmy,nu] = size(u);

  nf = 2*nu+1;
  uf = zeros(dmy,nf);
  uf(1) = (u(1)+u(2))/2;
  for i=2:nu-1
     uf(2*i)   = u(i);
     uf(2*i+1) = (u(i)+u(i+1))/2;
  end
  // strip old bc, add new bc
  uf = uf(2*nbc+1:nf-2*nbc);
  ufine = addbc(uf,nbc);
endfunction

// for a given source function we know the real solution; find rms error
function [e] = rmserror(u, x, nguard)
  nu = max(size(u));
  realsoln = (sin(4*x)/16);
  e = sqrt(sum((u(nguard+1:nu-nguard)-realsoln(nguard+1:nu-nguard))^2)/(nu-2*nguard));
endfunction
        
function [x,uo] = initialconditions(nx, h, nguard)
  n = nx + 2*nguard
  x = (1-nguard:1:nx+nguard)*h;
  
  uo = zeros(x);
endfunction


function [allt,x,allu,alle]=runredblack(nx,ndt)
  
  dx = 2.*%pi/nx;
  dt = dx*dx/(1.*16.);
  nguard = 2;
  
  [x,uo] = initialconditions(nx, dx, nguard);
  uo = setbc(uo,2);
  nu = max(size(uo));
  rhs = sourceterm(x);
  
  allt = 0:ndt;
  allu = zeros(ndt+1,nu);
  alle = zeros(ndt+1);
  
  allu(1,:) = uo
  alle(1) = rmserror(uo,x,nguard);
  for i=2:ndt+1 
    [du, tmp] = redblack(allu(i-1,:),nu,dx,dt,nguard,rhs,1.);
    allu(i,:) = setbc(tmp,nguard);
    alle(i)   = rmserror(allu(i,:),x,nguard);
  end
endfunction


function [allt,x,allu,alle]=runcoarsegrid(nx,ndt,nr)
  
  dx = 2.*%pi/nx;
  dt = dx*dx/(1.*16.);
  
  nguard = 2;
  
  [x, uo] = initialconditions(nx, dx, nguard);
  [dmy, nu] = size(uo);
  rhs = sourceterm(x);

  uo = setbc(uo,nguard);


  nsteps = ndt/(3*nr);
  ndt = nsteps*3*nr;
  
  allt = (0:1:ndt);
  allu = zeros(ndt,nu);
  
  allu(1,:) = uo;
  alle(1) = rmserror(uo,x,nguard);
  
  count = 2;

  for j=1:nsteps do
    
      // do nr relaxation sweeps at the finest level
      
      for i=1:nr do
          [resid, allu(count,:)] = redblack(allu(count-1,:),nu,dx,dt,nguard,rhs,1.);
          allu(count,:) = setbc(allu(count,:),nguard);
          allt(count) = allt(count-1)+1;
          alle(count) = rmserror(allu(count,:),x,nguard);
          count = count + 1;
      end

      // Now coarsen, including the residual
      [xcoarse,ucoarse] = initialconditions((nx+1)/2-1,2*dx,nguard)
      [dmy, nuc] = size(ucoarse);
      coarserhs = coarsen(resid,nguard);
      
      // and evolve the error equation on the coarse grid
      for i=1:2*nr do
        [resid, ucoarse] = redblack(ucoarse, nuc, 2*dx, 4*dt, nguard, coarserhs,1.);
        ucoarse = setbc(ucoarse,nguard);
      end
    
      // now update the fine grid solution
      correct = finer(ucoarse,nguard);
      allu(count,:) = allu(count-1,:) + correct;
      allt(count) = allt(count-1) + 2*nr;
      alle(count) = rmserror(allu(count,:),x,nguard);
      count = count + 1;
      
      // and relax another nr sweeps here.
      for i=1:nr do
          [resid, allu(count,:)] = redblack(allu(count-1,:),nu,dx,dt,nguard,rhs,2.);
          allu(count,:) = setbc(allu(count,:),nguard);
          allt(count) = allt(count-1)+1;
          alle(count) = rmserror(allu(count,:),x,nguard);
          count = count + 1;
      end
  end
  allt = allt(1:count-1);
  allu = allu(1:count-1,:);
  alle = alle(1:count-1,:);
endfunction