3738.sagews

import numpy as np

def f(t,y):
    return (-y**2)

def yip1(yi,h):
    delta=1+4*h*yi
    return((-1+sqrt(delta))/(2*h))

def Euler (a,b,N,y0):
         h=(b-a)/N
         y=np.empty(N+1)
         y[0]=y0
         for i in range (N) :
            t=a+i*h
            y[i+1]=yip1 (y[i],h) # Euler implicite
            # y[i+1]=y[i]+h*f(t,y[i]) # Euler explicite
         return (y)

m=Euler(0,5,50,1)

m

array([ 1. , 0.91607978, 0.84472393, 0.78335883, 0.73006006, 0.68336173, 0.64212879, 0.60546947, 0.57267392, 0.5431705 , 0.51649391, 0.49226175, 0.47015699, 0.44991467, 0.43131169, 0.41415893, 0.39829503, 0.38358155, 0.36989902, 0.35714385, 0.34522577, 0.33406577, 0.32359444, 0.3137505 , 0.30447971, 0.29573386, 0.28746996, 0.27964957, 0.27223821, 0.26520485, 0.25852151, 0.2521629 , 0.24610608, 0.24033022, 0.23481634, 0.22954716, 0.22450682, 0.21968086, 0.21505595, 0.21061988, 0.20636137, 0.20227006, 0.19833633, 0.19455131, 0.19090677, 0.18739508, 0.18400914, 0.18074236, 0.17758859, 0.17454209, 0.17159752])

list_plot(list(m))