% Convection - Diffusion - 2D; Dirichle BC on [0,1]*[0,1]
% compute a rhs vector for the central difference scheme
% 
function rhs=cd_rhs(h,nanr,flag,tita)
global cdeps

xcoor=h:h:(1-h);ycoor=h:h:(1-h);
% make them in a vector form
[x,y]=meshgrid([h:h:(1-h)],[h:h:(1-h)]);
nr2 = nanr^2;
Xcoor=reshape(x',nr2,1);
Ycoor=reshape(y',nr2,1);
u_x =cd_derx(Xcoor,Ycoor);
u_y =cd_dery(Xcoor,Ycoor);
u_xx=cd_derxx(Xcoor,Ycoor);
u_yy=cd_deryy(Xcoor,Ycoor);
B1 = cd_coefx(Xcoor,Ycoor,flag,tita);
B2 = cd_coefy(Xcoor,Ycoor,flag,tita);
rhs = -cdeps*u_xx-cdeps*u_yy+B1.*u_x+B2.*u_y;


% add the boundary conditions for the central differences case      
for j=1:nanr,
    for i=1:nanr,
        ind = (j-1)*nanr+i;
        xi=i*h;yj=j*h;
    if i-1<1,
       rhs(ind,1) = rhs(ind,1)+...
       (cdeps + cd_coefx(xi,yj,flag,tita)*h/2)/h^2*cd_sol_xy(xi-h,yj);
    end
    if i+1>nanr,
       rhs(ind,1) = rhs(ind,1)+...
       (cdeps - cd_coefx(xi,yj,flag,tita)*h/2)/h^2*cd_sol_xy(xi+h,yj);
    end
    if j-1<1,
       rhs(ind,1) = rhs(ind,1)+...
       (cdeps + cd_coefy(xi,yj,flag,tita)*h/2)/h^2*cd_sol_xy(xi,yj-h);
    end
    if j+1>nanr,
       rhs(ind,1) = rhs(ind,1)+...
       (cdeps - cd_coefy(xi,yj,flag,tita)*h/2)/h^2*cd_sol_xy(xi,yj+h);
    end
             
    end
end

return