// Definição do lado direito da EDO

m=1;g=10;l=1


function dz=f(t,z)
  dz(1)=z(2);
  dz(2)=-(g/l)*sin(z(1)); //-0.01*z(2)+0.1*sin(sqrt(10)*t);
endfunction

z_sol=RK4([%pi/4;0],0,80,f,40000);

plot(z_sol(1,1:100:$),z_sol(2,1:100:$));
xgrid;

o=z_sol(1,:);
vo=z_sol(2,:);

// A energia mecânica total deve ser uma constante.
// O vetor Emec guarda a energia mecânica total. 

Emec= (1/2)*m*(l*vo)^2 - g*l*cos(o);



scf();

// Podemos notar que a variação é pequena mas existe.

plot((Emec-Emec(1))/Emec(1));
xgrid;

