ComputeAx_collocatedVs2.m

% ComputeAx
function R=ComputeAx_collocatedVs2(X, rr, tt ,delr, delth, r_n, t_n,dt,Re)
% int ii r direction,jj theta direction; double P_rr, P_tt, P_r
%X=reshape(X1,r_n,t_n);

R=zeros(r_n,t_n);
for ii=2:r_n-1 %r direction
    for jj=1:t_n % theta direction
        % for theta derivative
        if(jj==1)
            
            Vss_tt = (X(ii,jj+1)+X(ii,t_n)-2*X(ii,jj))/(delth^2);
            
        elseif (jj==t_n)
            Vss_tt = (X(ii,1)+X(ii,jj-1)-2*X(ii,jj))/(delth^2);
        else
            Vss_tt = (X(ii,jj+1)+X(ii,jj-1)-2*X(ii,jj))/(delth^2);
        end
        
        % for first r derivatives
        if (ii==2)
            Vss_r=(X(ii+1,jj)-0)/(2*delr);% dirichlet
        elseif (ii==r_n-1)
            % At the farfield
            %              if (tt(ii,jj)<pi/2 || tt(ii,jj)>3*pi/2 ) % neumann
            %                  Vss_r=(X(ii,jj)-X(ii-1,jj))/(2*delr);
            %                  % downstream
            %             else
            %
            Vss_r=(0-X(ii-1,jj))/(2*delr);% Dirichlet
            %              end
            
        else
            Vss_r=(X(ii+1,jj)-X(ii-1,jj))/(2*delr); %
        end
        
        % for second r derivatives
        if (ii==2)
            Vss_rr=(X(ii+1,jj)-2*X(ii,jj)+0)/(delr^2);
            
        elseif (ii==r_n-1)
            
            %              if (tt(ii,jj)<pi/2 || tt(ii,jj)>3*pi/2 )
            %                 Vss_rr=(-1*X(ii,jj)+X(ii-1,jj))/(delr^2);
            %
            %            else
            %
            
            Vss_rr=(0-2*X(ii,jj)+X(ii-1,jj))/(delr^2);% Dirichlet
            %              end
            
        else
            Vss_rr=(X(ii+1,jj)-2*X(ii,jj)+X(ii-1,jj))/(delr^2);
        end
        
        
        R(ii,jj)= (dt/Re)*( Vss_rr +  (1/(rr(ii,jj)^2))*Vss_tt + (1/rr(ii,jj))*Vss_r ) ...
            -((dt/(Re*rr(ii,jj)^2)) +1)*X(ii,jj);
        
    end
end
%R(r_n-1,1:2)=X(r_n-1,1:2);
%R=R1(:);

end