function [dq]=dK(p),
    dq(1)=p(1)
endfunction
function [dp]=dV(q),
    dp(1)=q(1)
endfunction
		        
function [q,p]=pasEulerA(dK,dV,nv,h,qn,pn),
     dq=dK(pn);
     for v=1:nv,
        q(v)=qn(v)+h*dq(v)
     end;
     dp=dV(q);
     for v=1:nv,
        p(v)=pn(v)-h*dp
     end;
endfunction
		    
function [mq,mp,t]=eulerA(qi,pi,T,N,dK,dV,nv),
    h = T/N;
    mq = zeros(nv,N+1);
    mp = zeros(nv,N+1);
    t=zeros(1,N+1);
    mq(:,1)=qi;
    mp(:,1)=pi;
    t(1)=0;
    for k=1:N,
        [mq(:,k+1),mp(:,k+1)]=pasEulerA(dK,dV,nv,h,mq(:,k),mp(:,k));
        t(k+1)=t(k)+h;
    end;
endfunction
			
qi=[1];
pi=[0];
nv=1;
N=1000;
T=20;
[mq,mp,t]=eulerA(qi,pi,T,N,dK,dV,nv);
			
figA=scf();
plot2d(t,mq(1,:));
xtitle('','t','q');
			
T=20;
N=1000;
[mq,mp,t]=eulerA(qi,pi,T,N,dK,dV,nv);
		    
figB=scf();
plot2d(t,abs(mq(1,:)-cos(t)));
xtitle('','t','erreur');
			
T=2000;
N=T*50;
[mq,mp,t]=eulerA(qi,pi,T,N,dK,dV,nv);
			
figC=scf();
plot2d(mq(1,:),mp(1,:),rect=[-1.5,-1.5,1.5,1.5]);
xtitle('','p','q');
			
Em=zeros(1,N+1);
for k=1:N+1,
    Em(k)=0.5*mp(1,k)^2+0.5*mq(1,k)^2;
end;
			
figD=scf();
plot2d(t,Em,rect=[0,0,T,1]);
xtitle('','t','Energie');
			
