function yp=f(t,y)
// yp=f(t,y) désigne y'=f(t,y)
    v0=4.50
    g=9.81
    L=0.6
    c=(v0/L)^2-((2*g)/L)
    yp=sqrt(c+((2*g)/L)*cos(t))
endfunction

function y=Euler(f,a,b,y0,N)
// Données : f est la fonction définie par y'=f(t,y)
//          la fonction y étant définie sur l'intervalle [a,b]
//          y0 est la condition initiale y(a)
//          N est le nombre de points souhaités
// Résultat :   renvoie un vecteur ligne dont les composantes sont
//              les yk approchant y(a+kh) selon la méthode d'Euler
    y=[y0];
    t=a;
    h=(b-a)/N
    for i=1:N
        y=[y,y(i)+h*f(t,y(i))]
        t=t+h
    end
endfunction

a=0;
b=12;
y0=0;
N=200;
h=(b-a)/N;
temps=[a:h:b];
theta=[]
theta=Euler(f,a,b,y0,N);
thetap=[]
thetap=f(temps,theta)
clf
plot(theta,thetap,"r")
legend("Portrait de phase du pendule")
xlabel("position angulaire en rad")
ylabel("vitesse angulaire en rad/s")