import numpy as np, matplotlib.pyplot as plt 

# Resolver a equação diferencial u'(s)=u^2+2scos(u(s)), com u(0)=2.

def f(s,u):   # du/dt=2su
    p=u**2+2*s*np.cos(u)
    return p
    
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=0.5/n
u=[]
u.append(2)

t=np.linspace(0,0.5,n+1)

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