function draw()
for i=1:401
for j=1:401
x=-2+4*i/401;
y=-2+4*j/401;
s(i,j)=lyapunov2(x,y);
end
imagesc(s), drawnow % draw a map of lyapunov exponents
end
[x,y]=meshgrid(-2:0.01:2,-2:0.01:2);
X=2.*x./(1+x.^2+y.^2);
Y=2.*y./(1+x.^2+y.^2);
Z=(-1+x.^2+y.^2)./(1+x.^2+y.^2); % redraw on rieman sphere
figure
colormap hsv
[xx,yy,zz]=sphere(401);
props.EdgeColor = 'none';
props.Cdata = s;
surface(X,Y,Z,props);
drawnow , rotate3d on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function lambda = lyapunov2(x,y)
vec1=[1;0];
vec2=[0;1];
for i=1:50
x1=((y^3-1+(-1/5)*(x^3+y^3+1))*x)/((x^3-y^3+(17/15)*(x^3+y^3+1))); % a desboves map
y1=((1-x^3+(7/15)*(x^3+y^3+1))*y)/((x^3-y^3+(17/15)*(x^3+y^3+1)));
x=x1;
y=y1;
J=[-3/5*x^3/(32/15*x^3+2/15*y^3+17/15)+(4/5*y^3-6/5-1/5*x^3)/(32/15*x^3+2/15*y^3+17/15)-32/5*(4/5*y^3-6/5-1/5*x^3)*x^3/(32/15*x^3+2/15*y^3+17/15)^2 12/5*y^2*x/(32/15*x^3+2/15*y^3+17/15)-2/5*(4/5*y^3-6/5-1/5*x^3)*x/(32/15*x^3+2/15*y^3+17/15)^2*y^2;
-8/5*x^2*y/(32/15*x^3+2/15*y^3+17/15)-32/5*(22/15-8/15*x^3+7/15*y^3)*y/(32/15*x^3+2/15*y^3+17/15)^2*x^2 7/5*y^3/(32/15*x^3+2/15*y^3+17/15)+(22/15-8/15*x^3+7/15*y^3)/(32/15*x^3+2/15*y^3+17/15)-2/5*(22/15-8/15*x^3+7/15*y^3)*y^3/(32/15*x^3+2/15*y^3+17/15)^2];
% jacobian matrix
vec1=J*vec1;
vec2=J*vec2;
vec1 = vec1;
%vec2 = vec2 - (dot(vec1,vec2)/dot(vec1,vec1))*vec1; % GS
l1=norm(vec1);
%l2=norm(vec2);
h1(i)=log(l1)/i; % maximal lyapunov exponent
%h2(i)=log(l2)/i;
end
lambda=h1(end);
