Tmax=30;
h=10^(-2);
a=0.2;

Niter=Tmax/h;

O1=1; V1=0;
O2=0;   V2=0;

theta=zeros(4,Niter);
theta(:,1)=[O1;V1;O2;V2];

//-------Force-----------
function y=F(x)

y=zeros(4,1);
y(1)=x(2);y(3)=x(4);
y(2)=a*(x(3)-x(1))-sin(x(1));
y(4)=-a*(x(3)-x(1))-sin(x(3));

endfunction;
//-----------------------

for i=1:(Niter-1)

   //theta(:,i+1)=theta(:,i)+h*F(theta(:,i)); // EULER
   p1=F(theta(:,i));                    //RK4
   p2=F(theta(:,i)+h/2*p1);             //RK4
   p3=F(theta(:,i)+h/2*p2);             //RK4
   p4=F(theta(:,i)+h*p3);               //RK4  
   theta(:,i+1)=theta(:,i)+h*(1/6*p1+2/6*p2+2/6*p3+1/6*p4); //RK4
end;

xbasc(); x=0:h:(Tmax-h); 
plot2d([x;x]',theta([1,3],:)'); //trajectoire
// energie=1/2*(theta(2,:).^2+theta(4,:).^2)+2-cos(theta(1,:))-cos(theta(3,:))+1/2*a* ((theta(1,:)-theta(3,:)).^2) ; plot2d(x,energie,rect=[0,0,Tmax,1]);  //energie
