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,:);