function[R]=rhs_compact(Q,JM,KM,dx,dy,A,B,epse,jmax,kmax,M_inf) %B is the banded matrix B(1,4,1) used for the compact scheme R = zeros(size(Q)); [Qprime_x,Qprime_y] = diff_compact(Q,JM,KM,dx,dy,jmax,kmax,M_inf); for m = 1:3 nls R(JM,KM,m) = 0.0; for n = 1:3 R(JM,KM,m) = R(JM,KM,m) + ... A(m,n)*Qprime_x(JM,KM,n) + ... B(m,n)*Qprime_y(JM,KM,n) ; end end if epse ~= 0.0 Qx_d(JM,KM,:) = (Q(JM+1,KM,:) -2*Q(JM,KM,:) + Q(JM-1,KM,:))/dx; Qy_d(JM,KM,:) = (Q(JM,KM+1,:) -2*Q(JM,KM,:) + Q(JM,KM-1,:))/dy; R(JM,KM,:) = R(JM,KM,:)-epse*(Qx_d(JM,KM,:) + Qy_d(JM,KM,:)); end R(JM,KM,:) = -R(JM,KM,:);