%matplotlib inline
import matplotlib.pyplot as plt #Paquete para graficar
import numpy as np
from numpy import linalg
import scipy as sp
from scipy.optimize import broyden2

# condiciones iniciales
ya=1 # valor de la funcion en el borde
yb=3 # valor de la funcion en el borde
a=0 # inicio intervalo
b=1 # fin intervalo
n=100 # numero de intervalos utilizados en solver
h=(b-a)/n
solucion=np.zeros(n+1)

# BVP -> IVP
def Y_prima(u, y, x):
    Y = np.zeros(2)
    Y[0] = y**2+np.sin(x)
    Y[1] = u
    return Y

# Euler's method
def euler_method(a, b, n, Y0, Yprima):
    Yi = Y0
    for i in range(n+1):
        solucion[i]=Yi[1]
        Yi = Yi + h*Yprima(Yi[0], Yi[1], i*h+a)
    solucion[i]=Yi[1]
    return Yi


# Función error a la que debemos encontrar raíces
def F(s):
    Y_0=np.array((s, ya))
    yaprox=euler_method(a, b, n, Y_0, Y_prima)
    return yaprox[1]-yb


s=1
# Ejemplo de la funcion error
print("Error usando s={}: {}".format(s,F(s)))


# Root finder
print("Raíz encontrada por método Broyden: {}".format(broyden2(F,10)))
alpha=broyden2(F,10)

# Aplicación de solver con la raíz encontrada
Ysol=euler_method(a, b, n, np.array((alpha,1)), Y_prima)
print("Solución de BVP | y({})={} | y'({})={}".format(b,Ysol[1],b,Ysol[0]))

# Grafico solucion
plt.plot(np.arange(a,b+h,h), solucion,'r',label='p = 0')
plt.show()