function [x,y] = stablreg(a,b,c,n)
x=normrnd(0,1,n,1);
f = exp(10*x)./(1+exp(10*x));
y=a + b*f + c .* stabrnd(1.6, length(x));
y=y.^3;
plot(x,y,'.');
hold on
z=-2:.01:2;
f = exp(10*z)./(1+exp(10*z));
plot(z,(a+b*f).^3,'-');
hold off;
