function stability()
% Shows how the too large integration timestep may lead to
% numerical instability in a seemingly unconditionally stable system.
% Filip Piekniewski 2009
% works in Matlab and Octave
figure;
cm=jet;
i=0;
range=0.01:0.0005:0.0275;
range=wrev(range);
for tau=range;

i=ceil(((tau-min(range))*64)/(max(range)-min(range)));
if (i==0) i=1; end;
tspan=0:tau:0.5;

for v0=-4.4:1.3:10
v=[];
v(1)=v0;
for t=tspan
    v(end+1)=v(end)+tau*(-v(end)^3+3*v(end)^2-3*v(end)+1);
end;
plot(tspan,v(1:end-1),'color',cm(i,:));

hold on;
end;
axis([0 0.5 -5 10]);
drawnow;
end;

function y = wrev(x)
	if (nargin < 1|| nargin > 1); usage('wrev(x)'); end
	if(~isvector(x))
		error('x must be a vector');
	end
	l = length(x);
	k = 0:l-1;
	y = x(l-k);
