P-positions of Nim variant with multi-pile moves
Generally you compute the $N$ and $P$ positions for small numbers and hope to find a pattern. $(0,n)$ and (n,n) are $N$ because you take all the stones. $(1,2)$ is $N$ because you take all the stones. Keep working up to larger piles.
There can only be one $P$ position with a given number of stones in a pile because you can take stones from the other pile to get there. There can also only be one $P$ position with a given difference in the number of stones because the $(x,x)$ move will get you there.
I agree with your $(1,3),(2,6),$ and $(4,5)$ being $P$ positions. I think the next are $(7,10),(8,13),(9,16)$ but I may have missed something.
For an arbitrary impartial game like this (since the heaps aren't independent, I wouldn't really call this a heap game, but maybe a variation of Wythoff's game), I find it helpful to write some code to evaluate whether many positions are $\mathcal{N}$ or $\mathcal{P}$-positions (i.e. whether the next player to move has a winning strategy or not). Since this game has relatively simple rules, it can be coded up quickly, to see pictures of the $\mathcal{P}$-positions.
For example, the $\mathcal P$-positions for heap sizes up through 19 are shown here:
This shows that $(0,0)$, $(1,3)$, $(2,6)$, $(4,5)$, $(7,10)$, $(8,14)$, $(9,17)$, and $(13,18)$ are $\mathcal{P}$-positions.
No pattern is clear, so we can scale up to the heap sizes through 99:
Unfortunately, this doesn't seem to have a simple pattern. It's a little reminiscent of the pictures of positions of Wythoff's game whose Grundy value is equal to a fixed number. Those pictures have positions that are at most a bounded distance from the key lines of slope $\phi$ and $1/\phi$.
The best I can produce quickly is the heap sizes up to $399$:
They seem to be roughly hugging four lines, but there's some randomness - it's not as simple as the $\mathcal P$-positions of Wythoff's game.
Here is some Wolfram Language code (I promise no elegance or efficiency):
moves[pair_] :=
moves[pair] =
Union @@ Table[{If[pair[[1]] >= x, Sort@{pair[[1]] - x, pair[[2]]},
Nothing],
If[pair[[2]] >= x, Sort@{pair[[1]], pair[[2]] - x}, Nothing],
If[pair[[1]] >= x && pair[[2]] >= x,
Sort@{pair[[1]] - x, pair[[2]] - x}, Nothing],
If[pair[[1]] >= x && pair[[2]] >= 2 x,
Sort@{pair[[1]] - x, pair[[2]] - 2 x}, Nothing],
If[pair[[2]] >= x && pair[[1]] >= 2 x,
Sort@{pair[[2]] - x, pair[[1]] - 2 x}, Nothing]}, {x, 1,
Max @@ pair}];
pwins[pair_] := pwins[pair] = AllTrue[moves[pair], Not[pwins[#]] &];
MatrixPlot[
Table[If[pwins[Sort@{a, b}], 1, 0], {a, 0, 199}, {b, 0, 199}],
ColorFunction -> "Monochrome", DataReversed -> {True, False},
Frame -> False, ImageSize -> 800, PlotRangePadding -> 0]
You can run that code without Mathematica by going to the Wolfram Development Platform, clicking on "Create a New Notebook", pasting the code, and then using Shift+Enter or Numpad Enter to evaluate the code. The parameters in that code generate an image of the $\mathcal{P}$-positions for heap sizes up to 199.