lambda=-41
function xp=f(x)
    //lambda = -1
    xp = -x*x
endfunction

function dxp=df(x)
    //lambda=-1
    dxp = -2*x
endfunction

function x=sol(t)
    x=(1+t).^(-1)
endfunction

function x=newton_solve(y,dt)
    x=y
    for iter=1:10
        x=x-(x-dt*f(x)-y)/(1-dt*df(x))
        //disp(x-dt*f(x)-y)
    end
endfunction

function [tn,xn]=euler(n,x0)
    dt=1/n
    tn=linspace(0,1,n+1)
    xn=zeros(tn)
    xn(1)=x0
    for i=2:n+1
        xn(i)=xn(i-1)+dt*f(xn(i-1))
    end
endfunction

function [tn,xn]=euler2(n,x0)
    dt=1/n
    tn=linspace(0,1,n+1)
    xn=zeros(tn)
    xn(1)=x0
    for i=2:n+1
       xs=xn(i-1)+dt/2*f(xn(i-1))
       xn(i)=xn(i-1)+dt*f(xs)
    end
endfunction


function [tn,xn]=euler_imp(n,x0)
    dt=1/n
    tn=linspace(0,1,n+1)
    xn=zeros(tn)
    xn(1)=x0
    for i=2:n+1
        //xn(i)=xn(i-1)+dt*f(xn(i-1))
        xn(i)=newton_solve(xn(i-1),dt)
    end
endfunction



n=1000
clf
[tn,xn]=euler(n,1)
plot(tn,xn,color='green')
[tn,xn]=euler_imp(n,1)
errn=xn($)- sol(1)
[tn,xn]=euler2(n,1)
rk2_n=xn($)- sol(1)
plot(tn,xn,color='red')
plot(tn,sol(tn),color='blue')

//[tn2,xn2]=euler(2*n,1)
[tn2,xn2]=euler_imp(2*n,1)
err2n=xn2($)-sol(1)

[tn2,xn2]=euler2(2*n,1)
rk2_2n=xn2($)-sol(1)
disp(errn/err2n)
disp(rk2_n/rk2_2n)


