function [V_r,V_s] = GradVandermondeMatrices2D(N,r,s)

% Initialize the gradient of the modal basis at (r,s) at order N

N_p = length(r);
V_r = zeros(N_p); V_s = V_r;

[a,b] = rstoab(r,s);

col = 1;
for i = 0:N
	for j = 0:(N - i)
		[vect1,vect2] = GradModalBaseValues2D(a,b,i,j);
		V_r(:,col) = vect1'; V_s(:,col) = vect2';
		col = col + 1;
	end
end

end