Matplotlib error: Value Error: x and y must have same first dimension

[Abhishek]

I am trying to run some Python code for the last few hours. How can I achieve the effect of "dot divide" from Matlab, in the following code? I am having trouble working with list comprehension and numpy arrays and getting the following error:

    Traceback (most recent call last):
      File "Thurs.py", line 128, in <module>
        plt.plot(np.array(range(1,N/2+2)), Splot[alpha][iii,val]/utot[iii,val],color=cmap(iii/50))

    ValueError: x and y must have same first dimension

Code:

    import numpy as np
    import matplotlib as mpl
    import matplotlib.pyplot as plt
    from scipy.integrate import odeint
    
    N = 2
    K00 = np.logspace(3,5,101,10)
    len1 = len(K00)
    Qvec = np.logspace(-2,2,2,10)
    S10vec = np.logspace(2,6,2,10)
    len2 = len(Qvec)
    y0 = [0]*(3*N/2+3)
    Kplot = np.zeros((len1,len2))
    Pplot = np.zeros((len1,len2))
    S = [np.zeros((len1,len2)) for kkkk in range(N/2+1)]
    KS = [np.zeros((len1,len2)) for kkkk in range(N/2)]
    PS = [np.zeros((len1,len2)) for kkkk in range(N/2)]
    Splot = [np.zeros((len1,len2)) for kkkk in range(N/2+1)]
    KSplot = [np.zeros((len1,len2)) for kkkk in range(N/2)]
    PSplot = [np.zeros((len1,len2)) for kkkk in range(N/2)]
    
    for val in range(0,len2):
        for series in range(0,len1):
            K0 = K00[series]
            Q = Qvec[val]
            S10 = S10vec[val]
            r1 = 0.0001
            r2 = 0.001
            a = 0.001
            d = 0.001
            k = 0.999
            P0 = 1
            tfvec = [1e7, 1e10]
            tf = tfvec[val]
            time = np.linspace(0,tf,1001)
        
            def f(y,t):
                for alpha in range(0,(N/2+1)):
                    S[alpha] = y[alpha]
                for beta in range((N/2)+1,N+1):
                    KS[beta-N/2-1] = y[beta]
                for gamma in range(N+1,3*N/2+1):
                    PS[gamma-N-1] = y[gamma]
                K = y[3*N/2+1]
                P = y[3*N/2+2]
    
                ydot = np.zeros((3*N/2+3,1))
                B = range((N/2)+1,N+1)
                G = range(N+1,3*N/2+1)
                runsumPS = 0
                runsum1 = 0
                runsumKS = 0
                runsum2 = 0
    
                for m in range(0,N/2):
                    runsumPS = runsumPS + PS[m]
                    runsum1 = runsum1 + S[m+1]
                    runsumKS = runsumKS + KS[m]
                    runsum2 = runsum2 + S[m]
                    ydot[B[m]] = a*K*S[m]-(d+k+r1)*KS[m]
    
                for i in range(0,N/2-1):
                    ydot[G[i]] = a*P*S[i+1]-(d+k+r1)*PS[i]
    
                for p in range(1,N/2):
                    ydot[p] = -S[p]*(r1+a*K+a*P)+k*KS[p-1]+d*(PS[p-1]+KS[p])
    
                ydot[0] = Q-(r1+a*K)*S[0]+d*KS[0]+k*runsumPS
                ydot[N/2] = k*KS[N/2-1]-(r2+a*P)*S[N/2]+d*PS[N/2-1]
                ydot[G[N/2-1]] = a*P*S[N/2]-(d+k+r2)*PS[N/2-1]
                ydot[3*N/2+1] = (d+k+r1)*runsumKS-a*K*runsum2
                ydot[3*N/2+2] = (d+k+r1)*(runsumPS-PS[N/2-1])- \
                                a*P*runsum1+(d+k+r2)*PS[N/2-1]
    
                ydot_new = []
                for j in range(0,3*N/2+3):
                    ydot_new.extend(ydot[j])
                return ydot_new
            
            y0[0] = S10
            for i in range(1,3*N/2+1):
                y0[i] = 0
            y0[3*N/2+1] = K0
            y0[3*N/2+2] = P0
    
            soln = odeint(f,y0,time, mxstep = 5000)
            for alpha in range(0,(N/2+1)):
                S[alpha] = soln[:,alpha]
            for beta in range((N/2)+1,N+1):
                KS[beta-N/2-1] = soln[:,beta]
            for gamma in range(N+1,3*N/2+1):
                PS[gamma-N-1] = soln[:,gamma]
    
            for alpha in range(0,(N/2+1)):
                Splot[alpha][series,val] = soln[len(time)-1,alpha]
            for beta in range((N/2)+1,N+1):
                KSplot[beta-N/2-1][series,val] = soln[len(time)-1,beta]
            for gamma in range(N+1,3*N/2+1):
                PSplot[gamma-N-1][series,val] = soln[len(time)-1,gamma]
            
            u1 = 0
            u2 = 0
            u3 = 0
        
            for alpha in range(0,(N/2+1)):
                u1 = u1 + Splot[alpha]
            for beta in range((N/2)+1,N+1):
                u2 = u2 + KSplot[beta-N/2-1]
            for gamma in range(N+1,3*N/2+1):
                u3 = u3 + PSplot[gamma-N-1]
    
            K = soln[:,3*N/2+1]
            P = soln[:,3*N/2+2]
            Kplot[series] = soln[len1-1,3*N/2+1]
            Pplot[series] = soln[len1-1,3*N/2+2]
            utot = u1+u2+u3
    
        plt.figure(val)
        cmap = mpl.cm.autumn
        for iii in range(0,100,50):
            for alpha in range(0,(N/2+1)):
                plt.plot(np.array(range(1,N/2+2)), Splot[alpha][iii,val]/utot[iii,val],color=cmap(iii/50))
                plt.xlabel('i')
                plt.ylabel(r'$\frac{S_i}{S_{tot}}$ (nM)')
                plt.title('N = 20: Behavior at [S](0) = 10^' + str(log10(Qvec[val]) + 4) + '(nM)', fontsize=20)
        plt.show()

At the very least, can I extract the values that I need just for the plot?