for n = 1:2
    dz(2:end, 2:end) = dz(2:end, 2:end) + diff(hy(:, 2:JE)) - diff(hx(:, 2:JE));
    pulse = exp(-.5*(power((t0-T)/spread,2.0)));
    dz(ic,jc) = pulse;
    ez(1:IE,1:JE) = ga(1:IE,1:JE).*dz(1:IE,1:JE);
    hx(1:IE,1:JE-1) = hx(1:IE,1:JE-1) + 0.5*diff(ez, 1, 2);
    hy(1:IE-1,1:JE) = hy(1:IE-1,1:JE) + 0.5*diff(ez);
  end