I have a $2\times 3$ affine matrix $$ M = \pmatrix{a &b &c\\ d &e &f} $$ which transforms a point $(x,y)$ into $x' = a x + by + c, y' = d x + e y + f$
Is there a way to decompose such matrix into shear, rotation, translation,and scale ? I know there's something for $4\times 4$ matrixes, but for a $2\times 3$ matrix ?
You've written this somewhat unorthodoxly. To use that matrix for that transformation, one would more usually write
$$\pmatrix{x'\\y'\\1}=\pmatrix{a&b&c\\d&e&f\\0&0&1}\pmatrix{x\\y\\1}\;.$$
So the difference between a $2\times3$ matrix and a $4\times4$ matrix was only from your way of writing it; this works the same way as an affine transform in three dimensions, just with one fewer dimension. You can immediately factor out the translation,
$$\pmatrix{x'\\y'\\1}=\pmatrix{1&0&c\\0&1&f\\0&0&1}\pmatrix{a&b&0\\d&e&0\\0&0&1}\pmatrix{x\\y\\1}\;.$$
Then you just have to decompose $\pmatrix{a&b\\d&e}$ into shear, rotation and scaling in two dimensions.
[Edit in response to the comment:]
This isn't a unique decomposition, since you can do the shear, rotation and scaling in any order. Here's the decomposition I use:
$$A=\pmatrix{a&b\\d&e}=\pmatrix{p\\&r}\pmatrix{1\\q&1}\pmatrix{\cos\phi&\sin\phi\\-\sin\phi&\cos\phi}$$
with
$$ \begin{eqnarray} p&=&\sqrt{a^2+b^2}\;,\\ r&=&\frac{\det A}p=\frac{ae-bd}{\sqrt{a^2+b^2}}\;,\\ q&=&\frac{ad+be}{\det A}=\frac{ad+be}{ae-bd}\;,\\ \phi&=&\operatorname{atan}(b,a)\;, \end{eqnarray} $$
where $\operatorname{atan}$ is the two-argument arctangent function with operand order as in Java. This of course assumes $p\ne0$.