clear y1 y2;
y1(1) = 2.9;
y2(1) = 3.1;
nMax = 9;

for n=1:nMax
    y1(n+1) = 64 * y1(n) / (1 + y1(n))^3;
    y2(n+1) = 64 * y2(n) / (1 + y2(n))^3;
end;

n = 1:(nMax+1);
plot(n, y1, 'b-*', n, y2, 'r--x');
line([1, nMax+1], [3, 3], 'color', 'black', 'linewidth', 2, 'linestyle', ':');
title('y_{n+1}=64 y_n / (1 + y_n)^3: y_{eq}=3 is oscillatorily unstable');
xlabel('n'); ylabel('y_n');
legend('y_0 = 2.9', 'y_0 = 3.1', 'y_{eq} = 3', 'Location', 'Northwest');