Based on your example lines, you seem to have the further simplification that your delta x and delta y are uniform (in your example equal to 0.2). It's easy to generalize (and actually more readable) where dx~=dy
The finite difference coefficients matrix is very sparse, especially as N increases, so you want to use sparse matrices and assignments. Let's assume your coefficients matrix is called A rather than T_dot, with size [N,M].
First preallocate
A = sparse(N*M,N*M);
Then, kind of in the spirit of re-shaping, you can loop through the indices of A, and get the subscripts using ind2sub(). These subscripts allow you to easily find the elements that form your "cross", that is, in the internal nodes.
Once you find the nodes of interest using the subscripts, you can convert those back into the linear indices into the "unwrapped" version of A to do your assignments.
The most efficient way to populate sparse matrices (I think) is to vectorize the assignment. So you can create a temporary array s with the number of elements equal to the number of affected nodes in the "cross", i.e., 5, and write your equations into it.
Below example uses simple diffusion as an example, with D the diffusion coefficient. So rather than your central 1st order difference, it is based on 2nd order difference; should be easy to adapt to your case.
s = nan(1,5);
for k = 1:N*M
[i,j] = ind2sub([N,M],k);
% order of subscripts:
% [ i, j]
% [i+1, j]
% [i-1, j]
% [ i,j-1]
% [ i,j-2]
% this can only work for internal nodes so you need to have some if statements
% see below note on boundary conditions
% p will be linear the index into the grid
% which will correspond to the column index in the Jacobian
%
p = sub2ind([N,M],...
[i,i+1,i-1,i ,i ],...
[j,j ,j ,j-1,j-2]);
% p = [i,j]
s(1) = - D*2/dx/dx + D/dy/dy;
% p = [(i+1),j]
s(2) = + D/dx/dx;
% p = [(i-1),j]
s(3) = + D/dx/dx;
% p = [i,(j-1)]
s(4) = - D*2/dy/dy;
% p = [i,(j-2)]
s(5) = + D/dy/dy;
% subscript index into sparse A
A = A + sparse(k,p,s,N*M,N*M);
end
You didn't address boundary conditions in your question, but this scheme makes it pretty straightforward to deal with them as well. You can detect the corners and boundaries with if statements
for k = 1:N*M
[i,j] = ind2sub([N,M],k);
if i==1 && j==1 % corner
elseif i==1 % left boundary
elseif i==N % right boundary
elseif j==1 % top/bottom boundary
% ...etc
else % internal nodes
% then you can always do i-1 and i+1 and j-1 and j+1 without out-of-bounds errors
end
end
