function federpendel
T = 100; y_0 = 0; v_0 = 1;
s = 5; tau = 2^(-s); K = floor(T/tau);
y = zeros(K+1,2);
y(1,:) = [y_0,v_0];
for k = 1:K
    y(k+1,:) = y(k,:)+tau*Phi((k-1)*tau,y(k,:),tau);
    plot(tau*(0:k),y(1:k+1,1),'r');
    axis([0,T,-5,5]); drawnow;
end

function val = Phi(t,y,tau)
val = (f(t,y)+f(t+tau,y+tau*f(t,y)))/2;

function vec = f(t,y)
r = 1/10; D = 1; ell = 0;
vec = [y(2),-r*y(2)-D*(y(1)-ell)];