#!/usr/bin/python import numpy as np n=1000 i=np.linspace(0,n-1,n) [ij,ii]=np.meshgrid(i,i) A=10*np.eye(n)*(0.0001*(ij+ii)+1)+0.0001*(np.ones((n,n))-np.eye(n))*(ij+ii) f = 0.01*i x=np.linalg.solve(A,f)