from numpy import *
import numpy as np
from matplotlib import *
from scipy import *
from pylab import figure, show, setp
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

#We define a function which is going to be the recursive function.
def num_rossler(x1_n, y1_n, z1_n, x2_n, y2_n, z2_n, h, a, b, c, k, w1, w2):
    x1_n1=x1_n+h*(-w1*y1_n-z1_n+K*(x2_n-x1_n)) 
    y1_n1=y1_n+h*(w1*x1_n+a*y1_n)
    z1_n1=z1_n+h*(b+z1_n*(x1_n-c))   
 
    x2_n1=x2_n+h*(-w2*y2_n-z2_n+K*(x1_n-x2_n))
    y2_n1=y2_n+h*(w2*x2_n+a*y2_n)
    z2_n1=z2_n+h*(b+z2_n*(x2_n-c))
    return x1_n1, y1_n1, z1_n1, x2_n1, y2_n1, z2_n1
 
#Now we prepare some variables
#First the parameters
a=0.165
b=0.2
c=10
K=0.1
w1=0.99
w2=0.95

#Them the time interval and the step size
t_ini=0
t_fin=10000
h=0.01
numsteps=int((t_fin-t_ini)/h)
 
#using this parameters we build the time.
t=numpy.linspace(t_ini,t_fin,numsteps)
#And the vectors for the solutions
x1=numpy.zeros(numsteps)
y1=numpy.zeros(numsteps)
z1=numpy.zeros(numsteps)
x2=numpy.zeros(numsteps)
y2=numpy.zeros(numsteps)
z2=numpy.zeros(numsteps)
 
#We set the initial conditions
x1[0]=0.001
y1[0]=0.001
z1[0]=0.001
 
x2[0]=0.002
y2[0]=0.002
z2[0]=0.002
n=x1.size-1
 
#This is the main loop where we use the recursive system to obtain the solution
for k in range(0, n):
    #We use the previous point to generate the new point using the recursion
    [x1[k+1],y1[k+1],z1[k+1],x2[k+1],y2[k+1],z2[k+1]]=num_rossler(x1[k],y1[k],z1[k],x2[k],y2[k],z2[k],t[k+1]-t[k],a,b,c,K,w1,w2)