#-*-coding:utf-8-*-

from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import absolute_import

import matplotlib.pyplot as pl
import numpy as np
import math
from mpl_toolkits.mplot3d import Axes3D





def F(x,y,z):
	return 10*(y-x)
def G(x,y,z):
	return 28*x-y-x*z
def H(x,y,z):
	return x*y-(8/3)*z

def Euler2(F,G,H,a,b,x0,y0,z0,n):
	x=x0
	y=y0
	z=z0
	h=(b-a)/n
	Z=[z0]
	X=[x0]
	Y=[y0]
	while x<=b:
		xi,yi,zi=x,y,z
		x=xi+h*F(xi,yi,zi)
		y=yi+h*G(xi,yi,zi)
		z=zi+h*H(xi,yi,zi)
		Z.append(z)
		X.append(x)
		Y.append(y)
		x=xi+h
	z=z0
	y=y0
	x=y0
	while x>=a:
		xi,yi,zi=x,y,z
		x=xi-h*F(xi,yi,zi)
		y=yi-h*G(xi,yi,zi)
		z=zi-h*H(xi,yi,zi)
		Z.insert(0,z)
		X.insert(0,x)
		Y.insert(0,y)
		x=xi-h
	return X,Y,Z
	
a=0
b=50
n=10000
X,Y,Z=Euler2(F,G,H,a,b,0,0.1,0,n)
#pl.plot(Y,Z)
#pl.show()

fig = pl.figure()
ax = fig.gca(projection='3d')

ax.plot(X, Z, Y, label='orbite')
ax.legend()
pl.show()