For every boards are two chessboard patterns we can solve for. For example, for a 3x3 board we can have:

  101        010
  010   or   101
  101        010

where 1 is a white square and 0 is a black square.

We can represent any N×N board as a N² vector. So, we can represent one of the chessboard pattern as:

  010
  101 = 010101010 = y
  010

At the beginning of the game, the board is at some random initial pattern.

  110
  111 = 110111010 = y0
  010

The goal of the game is to find the necessary moves to form the chessboard pattern. Hence, we want to solve the equation

  (y - y0) = Mx

where M is a N²×N² matrix representing the effect of every possible moves. Pressing on a tile is equation to addition by 1 modulo 2.

  1 + 1 = 0 (mod 2)
  0 + 1 = 1 (mod 2)

For example, pressing the top left tile is the same as adding this effect vector to the board vector:

  110110000

If we do this for all 9 tiles of the 3x3 board, we find the matrix M to be:

  110110000
  111111000
  011011000
  110110110
  111111111 = M
  011011011
  000110110
  000111111
  000011011

where the i-th row is the effect of pressing the i-th tile. Now, we can just solve for x (assuming M is not singular). Re-using our above example, we find

  (y - y0) = Mx
  (010101010 - 110111010) = Mx
  100010000 = Mx
  111100100 = x   (by Gaussian elimination)

This means we can press the following tiles (in any order) on the board to solve it:

  111
  100
  100

If you are interested, I wrote some quick demo code that shows how this could be implemented in Python: https://gist.github.com/avassalotti/7452318c9e9b8dec637bc7c0...