function [R] = nl_rhs_cen(Q,JM,KM,dx,dy,A,B,epse)

%Get p,u, and v
p=Q(:,:,1);
u=Q(:,:,2)./Q(:,:,1);
v=Q(:,:,3)./Q(:,:,1);

sizeQ=size(Q);
N=sizeQ(1);
M=sizeQ(2);

A=zeros(N,M,3,3);
B=zeros(N,M,3,3);

A(:,:,1,1)=0;         A(:,:,1,2)=1;       A(:,:,1,3)=0;
A(:,:,2,1)=-u.^2+1;   A(:,:,2,2)=2*u;     A(:,:,2,3)=0;
A(:,:,3,1)=-u.*v;     A(:,:,3,2)=v;       A(:,:,3,3)=u;

B(:,:,1,1)=0;         B(:,:,1,2)=0;       B(:,:,1,3)=1;
B(:,
R = zeros(size(Q));
[Qx_c,Qy_c] = diff_cen(Q,JM,KM,dx,dy); 

for m = 1:3
  R(JM,KM,m) = 0.0;
  for n = 1:3
    R(JM,KM,m) = R(JM,KM,m) + ...
    A(m,n)*Qx_c(JM,KM,n) + ...
	B(m,n)*Qy_c(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,:);