from __future__ import division
from scipy import *
from scipy.integrate import odeint
from math import *

import pylab




x=arange(-0.5 , 1.5 , 0.01)
initial_y=[0.05] ## initial condition y(-0.5)=0.05
def f0(y,x):
    return y[0]**2+x**2

def f(y,x):
    y_dot=[0]
    y_dot[0]=f0(y,x)
    return y_dot
y=odeint(f,initial_y,x)


 #Plotting y vs x
pylab.figure(1)
pylab.plot(x, y[:, 0])
pylab.xlabel('x')
pylab.ylabel('y')
pylab.legend()
pylab.show()


# Direction field
# Creating a grid
rx = ravel(outer( ones(11) , arange(-0.5, 1.7, 0.2) ))
ry = ravel(outer( arange(-0.5, 1.7, 0.2), ones(11)  ))

x_component = zeros(len(rx))
y_component = zeros(len(rx))

# Caclulating the direction vector a each grid point, with rescaling factor a
a=1.0
for i in range(len(rx)) :
        x_component[i] = a*1.0
        y_component[i] = a*(ry[i]**2+ rx[i]**2)

# Plotting the vector field

pylab.xlabel('x')
pylab.ylabel('y')
pylab.quiver(rx, ry, x_component, y_component)
pylab.show()