Ffunction[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
Qsize=size(Q);
M=Qsize(1);
N=Qsize(2);

Bx=diag(4*ones(M-2,1))+diag(ones(M-3,1),1)+diag(ones(M-3,1),-1);
By=diag(4*ones(N-2,1))+diag(ones(N-3,1),1)+diag(ones(N-3,1),-1);
Qx_c=zeros(N,M,3);
Qprime_x=zeros(size(Q));Qprime_y=zeros(size(Q));
%To make the matrix multiplication make sense I take
for i=1:3
    Qt(:,:,i)=Q(:,:,i)';
    Qx_c(KM,JM,i)=0.5*(Qt(KM,JM+1,i)-Qt(KM,JM-1,i))/dx;
    Qprime_x(KM,JM,i)=6*inv(Bx)*Qx_c(KM,JM,i);
end

Bx=diag(4*ones(M-2,1))+diag(ones(M-3,1),1)+diag(ones(M-3,1),-1);
By=diag(4*ones(N-2,1))+diag(ones(N-3,1),1)+diag(ones(N-3,1),-1);
%Apply the boundary conditions to Q while maintaining zeros elsewhere
%[Qtempx,Qtempy]=bc8_compact(Q,JM,KM,jmax,kmax,M_inf);
Qprime_x=zeros(size(Q));Qprime_y=zeros(size(Q));

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(JM,KM,i)=6*inv(B)*Qx_c(JM,KM,i);
Qprime_y(JM,KM,i)=6*inv(B)*Qy_c(JM,KM,i);

end