
function a = diffus(dstab,a)

%  This subroutine applies a simple diffusion
%  operation to the array a. It leaves the endpoints
%  unchanged (assumes nothing about the
%  boundary conditions).

global nz

a(2:nz-1) = a(2:nz-1)+dstab*(a(1:nz-2)-2*a(2:nz-1)+a(3:nz));


% =-------------------------------------------------------------------------=


