function [laplace] = Laplacian(p,nz,nx,p_pre,alpha,vel)

laplace = zeros([nz nx]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%extrapolation
topOuter    = zeros([1 nx]);
leftOuter   = zeros([1 nz]);
bottomOuter = zeros([1 nx]);
rightOuter  = zeros([1 nz]);
for k = 2:nx-1
    topOuter(k)    = p(2,k) + (2.0/(alpha*vel(1, k)))*(p(1,k) - p_pre(1,k));%4*p(1,k) + 4*p(3,k) - p(4,k) - 6*p(2,k);
    bottomOuter(k) = p(nz-1,k) - (2.0/(alpha*vel(1, k)))*(p(nz,k) - p_pre(nz,k)); %4*p(nz,k) + 4*p(nz-2,k) - p(nz-3,k) - 6*p(nz-1,k);
end
for k = 2:nz-1
    leftOuter(k)  = p(k,2) + (2.0/(alpha*vel(k,1)))*(p(k,1) - p_pre(k,1));%4*p(k,1) + 4*p(k,3) - p(k,4) - 6*p(k,2);
    rightOuter(k) = p(k,nx-1) -(2.0/(alpha*vel(k,nx)))*(p(k,nx) - p_pre(k,nx)) ;%4*p(k,nx) + 4*p(k,nx-2) - p(k,nx-3) - 6*p(k,nx-1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 3:nx-2
    laplace(2,k) = 1/12.0*(16*p(2,k+1) - 30*p(2,k) + 16*p(2,k-1) - p(2,k-2) - p(2,k+2)) + ...
        1/12.0*(16*p(1,k) - 30*p(2,k) + 16*p(3,k) - p(4,k) - topOuter(k));
    
    laplace(nz-1,k) = 1/12.0*(16*p(nz-1,k+1) - 30*p(nz-1,k) + 16*p(nz-1,k-1) - p(nz-1,k-2) - p(nz-1,k+2)) + ...
        1/12.0*(16*p(nz,k) - 30*p(nz-1,k) + 16*p(nz-2,k) - p(nz-3,k) - bottomOuter(k));
end
for k = 3:nz-2
    laplace(k,2) = 1/12.0*(16*p(k,1) - 30*p(k,2) + 16*p(k,3) - p(k,4) - leftOuter(k)) + ...
        1/12.0*(16*p(k-1,2) - 30*p(k,2) + 16*p(k+1,2) - p(k-2,2) - p(k+2,2));
    
    laplace(k,nx-1) = 1/12.0*(16*p(k,nx) - 30*p(k,nx-1) + 16*p(k,nx-2) - p(k,nx-3) - rightOuter(k)) + ...
        1/12.0*(16*p(k-1,nx-1) - 30*p(k,nx-1) + 16*p(k+1,nx-1) - p(k+2,nx-1) - p(k-2,nx-1));
end

laplace(2,2) = 1/12.0*(16*p(2,1)-30*p(2,2)+16*p(2,3) -p(2,4)-leftOuter(2)) + ...
    1/12.0*(16*p(1,2)-30*p(2,2)+16*p(3,2) -p(4,2)-topOuter(2));

laplace(2,nx-1) = 1/12.0*(16*p(2,nx)-30*p(2,nx-1)+16*p(2,nx-2) -p(2,nx-3)-rightOuter(2)) + ...
    1/12.0*(16*p(1,nx-1)-30*p(2,nx-1)+16*p(3,nx-1)-p(4,nx-1)-topOuter(nx-1));

laplace(nz-1,2) = 1/12.0*(16*p(nz-1,1)-30*p(nz-1,2)+16*p(nz-1,3) -p(nz-1,4)-leftOuter(nz-1)) + ...
    1/12.0*(16*p(nz,2)-30*p(nz-1,2)+16*p(nz-2,2) -p(nz-3,2)-bottomOuter(2));

laplace(nz-1,nx-1) = 1/12.0*(16*p(nz-1,nx)-30*p(nz-1,nx-1)+16*p(nz-1,nx-2) -p(nz-1,nx-3)-rightOuter(nz-1)) + ...
    1/12.0*(16*p(nz,nx-1)-30*p(nz-1,nx-1)+16*p(nz-2,nx-1)-p(nz-3,nx-1)-bottomOuter(nx-1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 3:nz-2
    for j = 3:nx-2
        laplace(k,j) = 1/12.0*(16*p(k-1,j) - 30*p(k,j) + 16*p(k+1,j) - p(k-2,j) - p(k+2,j)) + ...
            1/12.0*(16*p(k, j-1) - 30*p(k,j) + 16*p(k, j+1) - p(k,j-2) - p(k,j+2));
    end
end
