clear;
h = 0.0001;
T = 16*pi;
N = ceil(T/h);

alpha = 1.0;

phi = zeros(N,1);
v = zeros(N,1);
t = zeros(N,1);

phi(1) = 1;
v(1) = 0;
t(1) = 0;

for j = 1:N-1
    
    t(j+1) = j*h;
    phi(j+1) = phi(j) + h*v(j);
    v(j+1) = v(j) + h* ( -alpha*v(j) -phi(j) );

end

plot(t,phi);
