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);
Qprime_x=zeros(size(Q));

M=Qsize(1);
N=Qsize(2);

Bx=diag(4*ones(N-2,1))+diag(ones(N-3,1),1)+diag(ones(N-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_xt=zeros(N,M,3);

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_xt(KM,JM,i)=6*inv(Bx)*Qx_c(KM,JM,i);
    Qprime_x(:,:,i)=Qprime_xt(:,:,i)';
end

Qprime_y=zeros(size(Q));
Qy_c=zeros(size(Q));
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
    
Qprime_y(JM,KM,i)=6*inv(By)*Qy_c(JM,KM,i);

end