How to solve systems of XOR equations?

Yes, you can use gaussian elimination to solve this. The key is to recognize that the XOR operation is equivalent to addition modulo 2. So the equation you wrote is equivalent to

i[0] = (p[0] + p[4] + ... ) mod 2

You can then set the whole system up as a matrix equation

M*p=i mod 2

You can solve this using Gaussian elimination as usual, except that all of your operations will be performed modulo 2. Since your matrix contains a lot of 0s, then you are going to have to use pivoting, but other than that, the algorithm is the same.