Cooley-Tukey Implementation of FFT in Matlab
The problem is in the definition of the twiddle factors. Your indices go from 1 to columns and from 1 to rows, whereas the indices in the power of the exponential function should go from 0 to columns-1 and from 0 to rows-1. So with this line of code in the twiddle factor loop everything should be fine:

matrix(y,x)=matrix(y,x)*exp((-2*pi*(y-1)*(x-1)/N)*1i);