If Statement Inside For Loop Creating Matrix

As stated in previous questions I am pretty new to python and I cannot figure out why my code is producing a matrix inside the for loop of my code. Well sort of, I think it has to do something with the if statements inside the for loop. I posted the code below. My code runs and produces almost the right output, but it creates a plot with many many lines created by the ArrayTheta and ArrayTheta1 arrays. Here is the picture: the correct image should look something like this:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

# Intialize the variables
g = 9.8
l = 9.8
step = .04
tmax = 60

theta = 0.2
theta1 = .001
ArrayTheta = [theta]
ArrayTheta1 = [theta1]

tx = np.arange(0,tmax,step)
omega = 0
omega1 = 0

q = 0.5
F_D = 1.2
w_d = float(2)/3
w_r = np.sqrt(g/l)

for i in np.arange(step,tmax,step):

    omega = omega - (g/l)*np.sin(theta)*step + -q*omega*step + F_D*np.sin(w_d*i)
    theta = theta + omega*step

    if theta < -np.pi:
        theta = theta + (2*np.pi)
    if theta > np.pi:
        theta = theta - 2*np.pi
    ArrayTheta.append(theta)


    omega1 = omega1 - (g/l)*np.sin(theta1)*step + -q*omega1*step + F_D*np.sin(w_d*i)
    theta1 = theta1 + omega1*step

    if theta1 > np.pi:
        theta1 = theta1 - 2*np.pi
    if theta1 < -np.pi:
        theta1 = theta1 + 2*np.pi
    ArrayTheta1.append(theta1)




plt.plot(tx,ArrayTheta)
plt.plot(tx,ArrayTheta1)
plt.show()

Now in the MatLab code a 1 x 1501 matrix is created and I was hoping to do the same in Python.

Edit

I included the MatLab code as well just so you can all see what is going on, and maybe it will be of use.

%This program creates two nonlinear driven pendulums with nearly identical 
%initial theta values and then plots the two. We then take the absolute value of the difference
%and plot this on a log plot. 

g=9.8; %Intialize gravity in m/s^2
l=9.8; %Intialize length of pendulum in meters
step=.04; %Intialize the step which is essentially change in time. 
tmax=60; %Intialize the tmax, which is how the program will run for in seconds. 
theta0=0.2; %Intialize theta_0 which we set to 0.2 radians. 

theta=theta0; %Intialize theta which will change over our programs duration. 
theta1=.001; %Intialize theta1 which will change over our programs duration. 
%To see that our generated data is correct we will create theta1=.001
%similiar to the textbook. 

ArrayTheta=theta; %Intialize the array which will contain the changing theta values.
ArrayTheta1=theta1;%Intialize the array which will contain the changing theta1 values.

tx=0:step:tmax; %tx runs over the duration of the program and will be used to help us plot later. 

omega0=0; %Intialize omega_0. Which is in radians/second.
omega=omega0; %Intialize omega which will change over the duration of the program. 
omega1=omega0; %Intialize omega1 which will change over the duration of the program.
q=.5; %strength of damping. 
%q1=2; %strength of damping.
F_D=1.2; %Initialize the driving force  
w_d=2/3; % Initialize the driving angular frequency
w_r=sqrt(g/l);%This helps us display resonance. 
           for t=step:step:tmax; 

        %The equation below this comment, computes omega over time with two
    omegaprev=omega; % We must include this step because we are using the Euler-Cromer method. 
    omega = omegaprev+((-g/l)*sin(theta)*step)+(-q.*omegaprev*step)+(F_D*sin(w_d*t))*step;
    theta = theta + (omega*step); % This function helps us find theta which is changing. 
    ArrayTheta= [ArrayTheta,theta]; %ArrayTheta stores all the changing values of theta through how many iterations the program runs. 


        %The equation below this comment, computes omega1 over time with two
    omegaprev1=omega1; % We must include this step because we are using the Euler-Cromer method. 
    omega1 = omegaprev1+((-g/l)*sin(theta1)*step)+(-q.*omegaprev1*step)+(F_D*sin(w_d*t))*step;
    theta1 = theta1 + (omega1*step); % This function helps us find theta1 which is changing. 
    ArrayTheta1= [ArrayTheta1,theta1]; %ArrayTheta1 stores all the changing values of theta through how many iterations the program runs. 

    %We need to create this while functions because we want the data to
    %stay in the range of -pi and pi. 
     while (theta < -pi); 

        theta = theta +(2*pi);             
     end
     while (theta > pi)
         theta = theta-(2*pi);      
     end
     %The set of while functions above made sure that the data for theta was staying
     %in the range of -pi and pi. 

      while (theta1 < -pi); 

        theta1 = theta1 +(2*pi);             
     end
     while (theta1 > pi)
         theta1 = theta1-(2*pi);      
     end
    %The set of while functions above made sure that the data for theta2 was staying
     %in the range of -pi and pi. 

     deltatheta=abs(ArrayTheta-ArrayTheta1); %We create this function to have a continuous log plot. 

           end

% We now want to plot our data. 

%First we will plot our theta values, theta & theta1 against one another. 
    figure
    hold on 
    plot(tx,ArrayTheta,'red');
    plot(tx,ArrayTheta1,'blue');
    xlabel('time (seconds)');
    ylabel('theta (radians)')
    hold off 

    %This plot, plots the log function to help us find the Lyapunov
    %constant.


    figure 
    semilogy(tx,deltatheta);
    xlabel('time (seconds)')
    ylabel('deltatheta')

    %We can say that the slope of the log function plot gives us 
    %the Lyapunov constant because log(deltatheta)= lambda * time. Where 
    %lambda is the Lyapunov constant. Using the basic fitting tool we get
    %that the slope is .0015 which means that the Lyapunov constant is
    %.0015.