Monte carlo works fine. Clearly you have a bug.
Why are you shuffling (with randperm) each suit before concatenating them? Would the order of the card in each suit matter when you draw one card at random?
Your allcards array is basically the numbers 1:13 repeated 4 times. You then draw one card (randsample is a bit overkill for this), so a number 1 to 13, and remove all instances of that number hence all 4 cards with the same value from the set, before putting back 3 of them at the end. Finally you draw one card again. Why such complicated logic? Why not draw two cards at once without replacement and be done with it?
A much simpler way to implement that simulation:
p = zeros(1, k);
allcards = repelem(1:13, 4);
for draw = 1:k
hand = allcards(randperm(52, 2));
p = diff(hand) == 0;
pp = cumsum(p) ./ (1:k);
which indeed converges to 100 x 3/51
edit: And I've finally spotted the bug in your code because when I looked at the p that your code generates it's all 0s except for the last element
if hands(1) == hands(2)
p(k) = 1;
with k constant and equal to the last element of p regardless of the loop index.