// Fonctions chute d'un corps


function Vprime=F(V);

Vprime = -k*V*V/m + g;


endfunction




function VV=eulerE(V0);
V=V0;
VV=V;
N=floor(T/dt);

for i=1:N,
    V=V+dt*F(V)
    VV=[VV,V];
end
endfunction

function VV=eulerI(V0 );
V=V0;
VV=V;
N=floor(T/dt);

for i=1:N,
    // on cherche le zero de la fonction 
    //  f(Y) = V(n) - dt *k*Y*Y/m + dt*g - Y
    //  V(n+1) = V(n) + dt F(V(n+1),m,k,g)
    // V(n+1) = V(n) - dt *k*V(n+1)*V(n+1)/m + dt*g
    // Y = V(n) - dt *k*Y*Y/m + dt*g
    
    deff('[y]=fsol1(x)','y=V - dt *k*x*x/m + dt*g - x');
    [xres]=fsolve([0.],fsol1);
    V=xres;
    VV=[VV,V];
end
endfunction

function VV=pointmilieu(V0 );
V=V0;
VV=V;
N=floor(T/dt);

for i=1:N,
    // on cherche le zero de la fonction 
     //  V(n+1) = V(n) + dt ( F(V(n+1)) + F(V(n))/2   
    // Y = V(n) - dt*k*Y*Y/m/2 + dt*g - dt*k*V(n)*V(n)/m/2 
     // Y + dt*k*Y*Y/m/2= V(n) + dt*g - dt*k*V(n)*V(n)/m/2 
    // -Y -dt*k*Y*Y/m/2 +V(n) +dt*g -dt*k*V(n)*V(n)/m/2 =0
    deff('[y]=fsol1(x)','y=V - 0.5*dt*k*V*V/m - 0.5*dt*k*x*x/m + dt*g - x');
    [xres]=fsolve([0.],fsol1);
    V=xres;
    VV=[VV,V];
end
endfunction



function V = solexacte()
   
    Vprime = sqrt(m*g/k);
     tprime = Vprime/g;
    time = [0:dt:dt*floor(T/dt)];
    V = Vprime*tanh(time/tprime);
endfunction

function Y=ff(V);

Y = -k*V*V/m + g;

endfunction


function [step,err] = erreurs()
    for i=1:5

        time = [0:dt:dt*floor(T/dt)];  
        VV1 = eulerE(V0);

        VV2 = eulerI(V0);

        VV3 = pointmilieu(V0);

        VV4 = solexacte();

        err(1,i) = abs(VV4(floor(T/dt))-VV1(floor(T/dt)));
        err(2,i) = abs(VV4(floor(T/dt))-VV2(floor(T/dt)));
        err(3,i) = abs(VV4(floor(T/dt))-VV3(floor(T/dt)));
        step(i) =   dt;
        dt = dt/2.;
    end

    plot2d(step,[err(1,:)',err(2,:)',err(3,:)'],[1,2,3],leg="EulerE@EulerI@point milieu)");
    
endfunction


