% pendulo simple sin rozamiento
set(0,'DefaultLineLineWidth',3);
clf;
% tiempo
t=[0:0.001:10];
% angulo inicial: probar pequenos (pi/6) y grandes (pi/1.1)
a0=pi/1.2;
% longitud
l0=1;
% velocidad angular inicial
v0=0;
% masa
m=2;
% coef. rozamiento: 0: no rozamiento (periodico)
r=0;
% gravedad
g=9.81;

% 2 Ley de Newton angular para el pendulo simple
% a es el angulo, w la velocidad angular
% a' = w
% w' = -g/l*sin(a) - r*l*w
f=@(t, X)...
    [X(2);...                % esto es a'=w
     -g/l0*sin(X(1)) - r*l0*X(2)];  % esto es v'=-g/l*sin(a)-r/m*w

% solucion
Y = heunvector(f,t,[a0;v0]);

% angulo en azul
% velocidad angular en rojo
plot(t, Y);
