Dear Friends, I have to solve a two dimensional coupled oscillator problem numerically. Hence the easiest way is to use the matrix representation of the creation and annihilation operator and then find the required hamiltonian. To find the time evolution of the state I have to multiply the state ψ with the unitary operator U=exp(-(i/h)HΔt). Where H is the coupled nonlinear harmonic oscillator hamiltonian. I did the problem with MATLAB using expm function. Getting better computational result I have to choose the size of the matrix H as 13225 by 13225. In matlab doing this integration for 1000 points will take more than 24 hours. Hence this is not a nice option. Hence I started programming using C gsl_libraries. Now the problem comes, how will I exponentiate a gsl_complex matrix ? that is the first difficulty. The second one is how will assign complex numbers to a gsl_matrix_complex. I am trying to diagonalize the hamiltonian then find the eigen values (D) and eigen vectors (V)then calculate V*Diag(i*dt*λ)*inv(V) . which will give me the expm result. But I dont know how to multiply a gsl_matrix with a complex number, do you have any ideas, regarding this problem ?