Stability of system of differential equations in Mathematica

Question

Can someone please run this Mathematica code on their computer/cloud with enough computation time available and post here the resulting image?

f[x_,y_]=-y*Sqrt[1-x^2]
g[x_,y_]=x*Sqrt[1-x^2]
i=0;
For[x0=-1.5,x0<=1.5,x0=x0+0.08,
For[y0=-1.5,y0<=1.5,y0=y0+0.08,
sol=NDSolve[{x'[t]==f[x[t],y[t]],y'[t]==g[x[t],y[t]],x[0]==x0,y[0]==y0},{x,y},{t,0,100}];
i=i+1;
t0=sol[[1,1,2,1,1,2]]*0.9;
While[Or[Abs[{x[t0]/.sol}[[1,1]]]>1000,Abs[{y[t0]/.sol}[[1,1]]]>1000],t0=t0/10];
a[i]=Show[ParametricPlot[{x[t],y[t]}/.sol,{t,0,t0},AxesOrigin->{0,0},PlotRange->{-1.5,1.5}],Graphics[{Arrowheads[0.03],Arrow[{{x[0],y[0]},{x[0.1],y[0.1]}}/.sol[[1]]]}]];
]]
Framed[Show[Table[a[m],{m,1,i}]]]

Note: If this question is deemed inappropriate for this site, please move it to the Mathematica stackexchange one.

Answer 1

Its not a good idea to use double for loops, even a single loop is very time consuming.

If you are interested in just phase portrait then use StreamPlot,

sp=StreamPlot[{-y*Sqrt[1 - x^2], x*Sqrt[1 - x^2]}, {x, -1, 1}, {y, -1, 1}, Axes -> True]

sol[y0_?NumericQ] := First@NDSolve[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]], 
                                   x[0] == y0, y[0] == y0}, {x, y}, {t, 0, 100}];

pp = ParametricPlot[Evaluate[{x[t], y[t]} /.sol[#] & /@ Range[-0.6, 0.6, 0.1]], {t,0,100}];

Show[pp, sp]