Problem with an algorithm to $3$-colour the edges of cubic graphs

Question

I'm currently trying to implement an algorithm to $3$-colour the edges of cubic graphs. (I want to do this with Matlab's Symbolic toolbox). After restricting to planar graphs to ensure the existence of a solution, I tried the following:

First, I observed that once a colouring is found the adjacency matrix $A$ split into three matrices with the following properties:

• $\hskip-5.0in A_1+A_2+A_3=A \tag{1}$
• $\hskip-5.0in A_k^2=1, \tag{2} $

where $1$ is the unit matrix and the idea is related to the concept of $1$-factors I think...

So I tried the following:

1. I set up matrices $X=(x_{kl})\odot A$, where $\odot$ is the Hadamard product ($Y$ and $Z$ accordingly),
2. plugged them into $(1)$ and $(2)$,
3. let the computer solve the set of equations
4. received some promising results for small dimensions (less than $10$)
5. and some frustating for larger ones (more than $20$).

I tried to use facts like, all entries in $A_k$ are either $0$ or $1$, so squares of certain variables in $(2)$ could be written as linear terms as well.

So my question is:

Is there a principle problem with this approach? Which other mathematical tricks can I use?