numpy - Python: Matrix multiplication error -

i found 2 matrices (mass , stiffness) using scipy.sparse.diags function , took 1 away other (mass - stiffness). when trying multiply new matrix vector u0 getting following

[ <6x6 sparse matrix of type '<type 'numpy.float64'>' 16 stored elements in compressed sparse row format>  <6x6 sparse matrix of type '<type 'numpy.float64'>' 16 stored elements in compressed sparse row format>  <6x6 sparse matrix of type '<type 'numpy.float64'>' 16 stored elements in compressed sparse row format>  <6x6 sparse matrix of type '<type 'numpy.float64'>' 16 stored elements in compressed sparse row format>  <6x6 sparse matrix of type '<type 'numpy.float64'>' 16 stored elements in compressed sparse row format>  <6x6 sparse matrix of type '<type 'numpy.float64'>' 16 stored elements in compressed sparse row format>] 

i don't understand result - getting 1d array of 6 matrices? if so, why getting , not simple vector?

the following code creating mass / stiffness matrix, finding u0 , dot product

import numpy np import math import scipy scipy.sparse import diags import scipy.sparse.linalg import matplotlib import matplotlib.pyplot plt   def mass_matrix(x0):     """finds mass matrix non uniform mesh x0"""     x0 = np.array(x0)     n = len(x0) - 1     h = x0[1:] - x0[:-1]     = np.zeros(n+1)     a[0] = 1     j in range(1,n):         a[j] = h[j-1]/3 + h[j]/3     a[n] = 1     b = h/6     b[0] = 0     c = h/6     c[-1] = 0     data = [a.tolist(), b.tolist(), c.tolist()]     positions = [0,1,-1]     mass_matrix = diags(data, positions, (n+1,n+1))      return mass_matrix  def poisson_stiffness(x0):     """finds poisson equation stiffness matrix nonuniform mesh x0"""     x0 = np.array(x0)     n = len(x0) - 1 # amount of elements; x0, x1, ..., xn     h = x0[1:] - x0[:-1]     = np.zeros(n+1)     a[0] = 1/h[0]     a[1:-1] = 1/h[1:] + 1/h[:-1]     a[-1] = 1/h[-1]     a[n] = 1/h[n-1]     b = -1/h     c = -1/h     data = [a.tolist(), b.tolist(), c.tolist()]     positions = [0, 1, -1]     stiffness_matrix = diags(data, positions, (n+1,n+1))      return stiffness_matrix  def nodal_quadrature_f(x0):     """finds nodal quadrature approximation of sin(pi x)"""     x0 = np.array(x0)     h = x0[1:] - x0[:-1]     n = len(x0) - 1     approx = np.zeros(len(x0))     in range(1,n):         approx[i] = math.sin(math.pi*x0[i])         approx[i] = (approx[i]*h[i-1] + approx[i]*h[i])/2      return approx  def initial_u(x0):      """finds initial u0"""     x0 = np.array(x0)     h = x0[1:] - x0[:-1]     n = len(x0) - 1     mass = mass_matrix(x0)     u0phi = nodal_quadrature_f(x0)     u0 = scipy.sparse.linalg.spsolve(mass, u0phi)      return u0  def timestepping(x0,theta,t,m):      x0 = np.array(x0)     h = x0[1:] - x0[:-1]     n = len(x0) - 1     dt = t/m     u0 = np.array(initial_u(x0))     mass = mass_matrix(x0)     stiffness = poisson_stiffness(x0)     iteration =  np.dot(mass - stiffness,u0)      print iteration 

i appreciate i've been stuck on few hours now.

edit: added information

the 6x6 matrix top recieved after call timestepping([0,0.2,0.5,0.6,0.9,1], 1, 1, 5)

i.e. x0 partition of doman 0 1. theta, t , m irrelevant havent added function yet.