import numpy as np, matplotlib.pyplot as plt 

# Resolver a equação diferencial u'(s)=(u1',u2'), com u(80)=(0,0). Estime u(10).

def f(s,u):   # 
    y=u[0]
    z=u[1]
    dy=(y-z)/s
    dz=-np.sqrt(1+dy**2)/2
    return np.array([dy, dz])
    
def phi(s,u):  # Fução de iteração do trapézio previsao-correcao-correcao
    p=u+h*f(s,u)
    v=u+h*( f(s,u)+f(s+h,p))/2
    w=u+h*( f(s,u)+f(s+h,v))/2
    return w

n=10
h=(10-80)/n
u=[]
y=[]
z=[]
u.append([0,0])
y.append(0)
z.append(0)

t=np.linspace(80,10,n+1)

for i in range(1,n+1):
    u.append(phi(t[i-1],u[i-1]))
    y.append(u[i][0])
    z.append(u[i][1])
    
    
print(f'Solucao aproximada u={u[n]}')
plt.plot(t,y,'o')
plt.plot(0*t,z,'o')
plt.show()