Thanks for all the effort...I think it will be useful...I will work it out.

well,this is just a part of the code..i'll be using these matrices..basically they are hamiltonian matrices on a grid...thanks anyway.

here is the code..i should have checked the preview..thanks for pointing out.

A:=proc(i, j,cc,numb)

    global mat1r,mat1i:
    mat1r[i,j]  := mat1r[i,j]  + coef*ec(i,j) - coef*((1/delz2) + (1/dely2) + (1/delx2)):       
    if j+8<(numbr*numbr*numbr*8) then
    mat1r[i,j+8] := mat1r[i,j+8] + 2*numb*coef/delz2
    fi:
    if (j+(numbr*8))<(numbr*numbr*numbr*8) then
    mat1r[i,j+(numbr*8)] := mat1r[i,j+(numbr*8)] + 2*numb*coef/dely2
    fi:
    if (j+(numbr*numbr*8))<(numbr*numbr*numbr*8) then
    mat1r[i,j+(numbr*numbr*8)] := mat1r[i,j+(numbr*numbr*8)] + 2*numb*coef/delx2
    fi:
    if (j+16)<(numbr*numbr*numbr*8)  then
    mat1r[i,j+16] := mat1r[i,j+16] - numb*coef/delz2
    fi:
    if (j+(numbr*16))<(numbr*numbr*numbr*8)  then
    mat1r[i,j+(16*numbr)] := mat1r[i,j+(16*numbr)] - numb*coef/dely2
    fi:
    if (j+(numbr*16*numbr))<(numbr*numbr*numbr*8)  then
    mat1r[i,j+(16*numbr*numbr)] := mat1r[i,j+(16*numbr*numbr)] - numb*coef/delx2
    fi:

end:
mat1r:=Matrix(8000000,8000000,storage=sparse,datatype=float[8],order=C_order):
mat1i:=Matrix(8000000,8000000,storage=sparse,datatype=float[8],order=C_order):

delx:=0.0015:dely:=0.0015:delz:=0.0015:
    numbr1:=100:chg:=0.3029:coef:=1e-6:
    numbr := numbr1:
    ncol:=numbr*numbr*numbr*8:
    delx2 := delx*delx:
    dely2 := dely*dely:
    delz2 := delz*delz:
    n3:=numbr*numbr*numbr:
    seq(A(0+(8*i),1+(8*i),0,1),i=1..n3):