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+Qtempx; Qy_c=Qy_c+Qtempy; Qprime_x(:,:,i)=6*inv(B)*Qx_c(:,:,i); Qprime_y(:,:,i)=6*inv(B)*Qy_c(:,:,i); end