Dear Numerical Python users, A small mistake when I cut and pasted the program: Here is the program again static PyObject* matrix(PyObject *self, PyObject *args){ double x,y; int size; int M,n,m; PyArrayObject *result; int dimensions[2]; Py_complex *data; Py_complex p; if(!PyArg_ParseTuple(args,"idd",&M,&x,&y)) return NULL; dimensions[0] = M; dimensions[1] = M; data = calloc(M*M+1,sizeof(Py_complex)); for(n=0;n < M;n++){ for(m=0; m<M; m++){ p.real=x; p.imag=y; data[n*M+m] = p; } } result = (PyArrayObject*)PyArray_FromDims(2,dimensions,PyArray_CDOUBLE); memcpy(result->data,data,M*M*sizeof(Py_complex)); return PyArray_Return(result); free(data); }