title: Speeding up Python when using nested loops
body: <p>This code, which contains 3 nested loops, in Matlab takes 30 sec to get run, instead, in Python more than 5 minutes. Any suggestion on how to make Python faster is highly appreciated! The code creates response spectra by using Newmark numerical solution. 
  I am aware that vectorizing will resolve the challenge but I could not find space for implementing it.  </p>
  
  <pre><code>from numba import jit
  import numpy as np
  import matplotlib.pyplot as plt
  import scipy as sp
  import timeit
  
  start = timeit.default_timer()
  sig=np.loadtxt('LXR1N.txt')
  t=sig[:,0]
  ag=0.01*sig[:,1]
  h=t[1]-t[0]
  f=-ag
  alfa=1/6; delta=0.5
  alfa0=1/(alfa*h**2)
  alfa1=delta/(alfa*h)
  alfa2=1/(alfa*h)
  alfa3=1/(2*alfa)-1
  alfa4=delta/alfa-1
  alfa5= h/2*(delta/alfa-2)
  alfa6= h*(1-delta)
  alfa7= delta*h;
  D=np.array([0.01,0.04,0.06,0.1])
  T=np.zeros(1198,float)
  d=np.zeros([len(t)])
  len(d)
  v=np.zeros([len(t)])
  a=np.zeros([len(t)])
  fb=np.zeros([len(t)-1])
  sd=np.empty((len(D),len(T)))
  sa=np.empty((len(D),len(T)))
  sv=np.empty((len(D),len(T)))
  aex=np.empty((len(D),len(T)))
  vex=np.empty((len(D),len(T)))
  for i in range(0,len(T)):
          T[i]=0.015+0.005*i
  i=None
  s=-1
  for q in D:
      s+=1
      for i in range(0,len(T)):
          omega=2*np.pi/T[i]
          c=2*q*omega
          k=omega**2
          kb=k+alfa0+alfa1*c
          d[0]=0
          v[0]=0
          a[0]=0
          for j in range(0,len(t)-1):
              fb[j]=f[j]+alfa0*d[j]+alfa2*v[j]+alfa3*a[j]+c*alfa1*d[j]+c*alfa4*v[j]+c*alfa5*a[j] 
              d[j+1]=fb[j]/kb
              a[j+1]=alfa0*d[j+1]-alfa0*d[j]-alfa2*v[j]-alfa3*a[j]
              v[j+1]=v[j]+alfa6*a[j]+alfa7*a[j+1]
          atot=a+ag
          atot1=-omega**2*d-2*q*omega*v
          sd[s,i]=max(abs(d))   
          sa[s,i]=omega**2*sd[s,i]
          sv[s,i]=omega*sd[s,i]
          aex[s,i]=max(abs(atot))
          vex[s,i]=max(abs(v))
  del i
  del j
  plt.figure(1,figsize=(10, 4))
  for i in range(sd.shape[0]):
      plt.plot(T,sd[i,:],label=str(D[i]*100)+"%")
      plt.xlabel("Period [s]")
      plt.ylabel("Sd [m]")
  plt.xlim([0,6])
  plt.ylim([0,0.6])
  plt.legend()
  plt.savefig('Sd.jpg')
  plt.show()
  plt.figure(2,figsize=(6, 3))
  for j in range(sa.shape[0]):
      plt.plot(T,1/9.81*sa[j,:],label=str(D[j]*100)+"%")
      plt.xlabel("Period [s]")
      plt.ylabel("Sa [g]")
  plt.xlim([0,6])
  plt.ylim([0,2.5])
  plt.legend()
  plt.savefig('Sa.png')
  plt.show()
  stop = timeit.default_timer()
  execution_time = stop - start
  
  print("Program Executed in: ",execution_time, 'seconds')
  </code></pre>