%%%%% Begin file ddz2.m
function [d2udz2] = ddz2(u)                 % DG second derivative operator
Globals;
du = zeros(2*K,1); du(:) = (u(vmapM)-u(vmapP))/2.0;
uin  = -u(1); du(1) = (u(1) -  uin )/2.0;
q = Dz*u - LIFT*(Fscale.*(nz.*du));         % auxiliary variable q = du/dz
dq = zeros(2*K,1); dq(:) = q(vmapM)-q(vmapP);
qout = -q(end);   dq(end) = (q(end)-qout);  % Neumann condition at outflow
tau = 1.0; fluxq = nz.*(dq/2.0+tau*nz.*du); % Stabilized central flux
d2udz2 = Dz*q - LIFT*(Fscale.*fluxq);
end
%%%%% End file ddz2.m