If you want it fast, you will have to define the matrix once by hand. Number the unknowns
U1 = u_{0,0}, U2 = u_{0,1},...,U(N+1)^2 = u_{N,N}
Proceed the same way with the equations EQN1,EQN2,...,EQN(N+1)^2 for U1,U2,...,U(N+1)^2.
Then your matrix will have an entry C in position (i,j) if the coefficient in front of Uj in EQNi equals C.
A slower, but more comfortable way is to write your equations from above as a function:
function res = fun(u)
and define
res(nx,ny) = E*u(nx,ny) - (-i*c*sqrt((nx+1)*ny)*u(nx+1,ny-1)+i*c*sqrt(nx*(ny+1))*u(nx-1,ny+1))
Then use "jacobianest" from the file exchange (or some faster routine) to calculate the Jacobian matrix A of the
function defined in "fun".
