Aqui está um exemplo de um programa em Octave que implementa o método de Heun para integrar a equação diferencial que mencionamos anteriormente e cria um gráfico dos resultados usando a biblioteca Gnuplot:

% Constantes
m = 2.6;  % Massa em g
L = 1.0;  % Comprimento em m
R = 0.03; % Raio em m

% Valores iniciais
theta0 = 0.05;  % Ângulo inicial em rad
w0 = 0.0;       % Velocidade angular inicial em rad/s
t0 = 0.0;       % Tempo inicial em s
h = 0.1;        % Passo de tempo em s

% Listas para armazenar os resultados
t = t0;
theta = theta0;
w = w0;

% Loop de tempo
while (t(end) < 100.0)
    % Cálculo dos valores intermediários
    k1x = h * w(end);
    k1v = h * (-1.0 if w(end) < 0.0 else 1.0);
    k2x = h * (w(end) + k1v);
    k2v = h * (-1.0 if w(end) + k1v < 0.0 else 1.0);

    % Atualização dos valores de θ e w
    theta = [theta theta(end) + (k1x + k2x) / 2.0];
    w = [w w(end) + (k1v + k2v) / 2.0];

    % Atualização do tempo
    t = [t t(end) + h];
end

% Cria o gráfico usando o Gnuplot
plot(t, theta, "r", "linewidth", 2);
hold on
plot(t, w, "b", "linewidth", 2);
legend("θ", "w");
title("Heun method");
xlabel("t (s)");
ylabel("θ (rad), w (m/s)");
grid on;