#
#octave_mechanics_singlependulum.m
#https://ix86.org/octave/index-main.html?octave_mechanics_singlependulum
#
function out = octave_mechanics_singlependulum(m,l,g,x0,x1)
dt =  0.01;
%m = 1;
%l = 0.5;
%g = 9.8;
%x_k = [30/180*pi ;
%          0        ];
x_k = [x0 ;
       x1];
x = [];
tim =[];

for t = 0 : dt : 10
    k1 = [ x_k(2);
          -g/l * sin(x_k(1)) ] * dt;
    x_tmp = x_k + k1/2;
    k2 = [ x_tmp(2);
          -g/l * sin(x_tmp(1)) ] * dt;
    x_tmp = x_k + k2/2;
    k3 = [ x_tmp(2);
          -g/l * sin(x_tmp(1)) ] * dt;
    x_tmp = x_k + k3;
    k4 = [ x_tmp(2);
          -g/l * sin(x_tmp(1)) ] * dt;
    x_k1 = x_k + (k1 + 2 * k2 + 2 * k3 + k4 ) /6 ;
    x = [x x_k];
    tim = [tim t];
    x_k = x_k1;
endfor

plot(tim, x(1, :),tim, x(2, :));
out=0;
endfunction

