function dy = Turing_ode(t,y,M,D,N,dx)

dy = 0*y;

A = y(1:N);
B = y((N+1):end);

diffmatrix = toeplitz([-2 1 zeros(1, N-2)]);
diffmatrix(end,1) = 1;
diffmatrix(1,end) = 1;


dy(1    :  N) = M(1,1) * A + M(1,2) * B + D(1) * diffmatrix * A / dx^2 - 10^(-5) * A.^3;
dy((N+1):end) = M(2,1) * A + M(2,2) * B + D(2) * diffmatrix * B / dx^2 - 10^(-5) * B.^3;