When I step through the code in ode45, x is a 7200 x 1, 1-d array. How do I ensure x remains a 2-d matrix?
To do that, write your own version of ode45. The Mathworks version will always pass in the boundary conditions as a column vector.
Most people deal with the situation by simply reshape() the input almost immediately.
Also what are the dimensions of y, when the solver has finished?
y will be length(t1) by numel(x).
Your function is required to return a column vector that is numel(x) by 1, and integrated versions of that column will become rows of the output.
What is it in the case x is a 2-d matrix
Then y will be length(t1) by numel(x).
You can reshape y afterwards, such as
yperm = permute(reshape(y, length(t1), size(x,1), size(x,2)),[2 3 1]);
and now yperm(:,:,iteration) will be the 2D results at time t1(iteration)