3736.sagews

import numpy as np

def yr(z,yi):
    delta=1+4*z*yi
    return ((-1+sqrt(delta)))

#Euler Implicite

def E_I(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]=yr(h,y[i])
    return (y)

E_I(0,5,5,1)

array([ 1. , 3.58257569, 7.52358574, 11.30738457, 14.07141968, 15.80560602])