
function z=f(t,y)
    z=-y/(1+t^2)
endfunction

function R=EulerExp(f,y0,t0,tmax,dt)
    T=[0]
    Y=[y0]
    for i=1:(tmax-t0)/dt
        T=[T,i*dt]
        Y=[Y,Y($)+dt*f(T($-1),Y($))]
    end
    R=[T;Y]
endfunction

Re=EulerExp(f,5,0,10,0.1)

scf(0)
clf()
plot2d(Re(1,:),Re(2,:),4)


function R=EulerImp(y0,t0,tmax,dt)
    T=[0]
    Y=[y0]
    for i=1:(tmax-t0)/dt
        T=[T,i*dt]
        Y=[Y,Y($)/(1+dt/(1+T($)^2))]
    end
    R=[T;Y]
endfunction

Ri=EulerImp(5,0,10,0.1)

plot2d(Ri(1,:),Ri(2,:),3)

function y=yexact(t)
    y=5*exp(-atan(t))
endfunction

plot2d([0:0.1:10],yexact([0:0.1:10]))

scf(1)
clf()

plot2d([0:0.1:10],abs(Re(2,:)-yexact([0:0.1:10])),3) //  courbe verte
plot2d([0:0.1:10],abs(Ri(2,:)-yexact([0:0.1:10])),4) // courbe bleue


