dx1dt(q,1) = k(2)*(1/(1+(x2/K(2))^n)) + k(1)*((x1/K(1))^n/(1+(x1/K(1))^n)) - b(1)*x1;
x2 is a 31 x 31 matrix, so (1+(x2/K(2))^n) is a 31 x 31 matrix. You have 1 / that. The / operator is "matrix right divide", formal name mrdivide(), and A/B is approximately equivalent to A * pinv(B) but with a requirement that the number of columns in A and B must be the same (the number of rows does not need to be the same.) The left side, 1, has one column, but the right side is 31 x 31, which is not the same number of columns.
You should consider using ./ instead of / if you mean element-by-element division.
Later in the expression you have k(1)*((x1/K(1))^n/(1+(x1/K(1))^n)) . The first part, k(1)*((x1/K(1))^n, is 31 x 31. The second part, 1+(x1/K(1))^n, is also 31 x 31. Because those have the same number of columns, the / operator is valid, and will act much like left side matrix-multiply by the inverse of the right side, doing a series of linear fits. The result would be 31 x 31.
However, if you use the / operator in that meaning, it looks suspicious compared to using ./ in the first part in a construct that is very similar. It is possible that is what you mean, but it is questionable, and leads one to ask whether perhaps this second one should be ./ as well, going for element-by-element division.
The third part of the expression is -b(1)*x1 which would be 31 x 31 . You would thus be adding or subtracting three 31 x 31 matrices, and that would be meaningful, giving a 31 x 31 result.
Now on the left hand side of the assignment, you have dx1dt(q,1) . However, q is undefined.
If we hypothesize that q is a scalar, then you have the problem that the left hand side would designate a single scalar location but the right hand side is 31 x 31, and 31 x 31 will not fit in a scalar location.
If we hypothesize that q is a vector or matrix that just happens to have (31*31) = 961 elements in it, then even though the left side and right side would both have 961 elements, with you using two indices, dx1dt(q,1) rather than one index dx1dt(q), then MATLAB will notice that the number of rows and columns do not match and will not permit the assignment.
If q does just happen to be a vector or array of 961 index values, then you would have to reshape() the right hand side into a vector in order for the assignment to work.
In such a case, dx1dt would have to be a column vector with at least 961 elements, or dx1dt would have to be a 2D array with at least 961 rows and an unknown number of columns.
In either situation, you would get a mismatch in sizes when you tried to plot the 31 x 31 x1 and x2 against the 961 or more rows by some number of columns dx1dt .
I would suggest to you that you do not want the (q,1) indexing in the assignment to dx1dt .