function[Qprime_x,Qprime_y]=diff_compact(Q,JM,KM,dx,dy,jmax,kmax,M_inf)

%B is the banded matrix B(1,4,1) used for the compact scheme
B=zeros(size(Q));
Qsize=size(Q);
M=Qsize(1);
N=Qsize(2);
B=diag(4*ones(M,1))+diag(ones(M-1,1),1)+diag(ones(M-1,1),-1);

%Apply the boundary conditions to Q while maintaining zeros elsewhere
%[Qtempx,Qtempy]=bc8_compact(Q,JM,KM,jmax,kmax,M_inf);

Qx_c=zeros(size(Q)); Qy_c=zeros(size(Q));
Qx_c(JM,KM,:) = 0.5*(Q(JM+1,KM,:) - Q(JM-1,KM,:))/dx;
Qy_c(JM,KM,:) = 0.5*(Q(JM,KM+1,:) - Q(JM,KM-1,:))/dy;

%Qx_c=Qx_c+Qtempx;
%Qy_c=Qy_c+Qtempy;

for i=1:3
    
%Qx_c(:,:,i)=Qx_c(:,:,i)+Qtempx(:,:,i);
%Qy_c(:,:,i)=Qy_c(:,:,i)+Qtempy(:,:,i);
Qprime_x(:,:,i)=6*inv(B)*Qx_c(:,:,i);
Qprime_y(:,:,i)=6*inv(B)*Qy_c(:,:,i);

end