// Sistema pêndulo-mola

// Constantes do sistema
g=10;m=1;k=30;L=1

// Definição do lado direito da EDO

function dz=f(t,z)
  dz(1)=z(2)
  dz(2)=-g/z(3)*sin(z(1))-2*z(2)/z(3)*z(4)
  dz(3)=z(4)
  dz(4)=g*cos(z(1 ))-k/m*(z(3)-L)+z(3)*z(2)^2
endfunction

// Condição Inicial
z0=[%pi/4;0;1;0]  // Pêndulo em repuso a 45º, mola não tensionada.

// Integração numérica da EDO no intervalo t, (0,20) com intervalo de tempo
// 0,001.
h=0.001
n=(20-0)/h
sol=zeros(4,n+1)
sol(:,1)=z0
t=0
for i=1:n
    sol(:,i+1)=RK4(sol(:,i),t,t+h,f)
    t=t+h
end

// Variáveis introduzidas para guardar o resultado da integração
 o=sol(1,:)
vo=sol(2,:)
 l=sol(3,:)
vl=sol(4,:)
clear sol  // como o resultado está guardado nas variáveis o,vo,l e vl, podemos 
           // liberar memória limpando a variável sol.

// A energia mecânica total é constante (no caso da solução exata).
// O vetor Emec guarda a energia mecânica. 

Emec=m/2*((l.*vo)^2+vl^2) -m*g*l.*cos(o) + 1/2*k*(l-L)^2

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

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

