As far as I understand the fist part of your question is regarding the options of the ode15s, and which ones you can use to make the solver more efficient/accurate.
Providing the Jacobian to the solver, especially if its reasonably easy to calculate, is a good idea to make things faster and more accurate.
Using sparse systems for small systems of equations ...
Some things I can think of:
use sparse matrices for Matrix1 and Matrix2 to speed up the computations of dZand dY
use larger integration tolerances reltoland abstol, especially if your are searching for steady-state solutions and/or do not need a precise resolution of the transient dynamics of your system.
you are solving with ode15s, which is an implicit ...
This post is better suited for Stackoverflow I think.
Anyway you can simply solve your problem by changing psi[:,ig]to psi[:,ig:ig+1]. Then the left-hand side is truly a nx1 matrix, not just a vector of size n.
Or you can remove the line nodes = np.array([nodes]).T which is useless here, and causes Numpy to transform the array "nodes" (size n), ...