Hi,

I want to run this code, but I do not draw diagrams. You can help me?

>restart;`read`(orbit.sav)

> dsnumsort := proc (numpr::list, Coor::list)

local i, j, n; global C1, C2, C3, V1, V2, V3;

n := nops(Coor); print("Order of the variables:");

for i from 2 to 2*n+1 do

for j to n do

if [numpr[i]] = select(has, numpr, diff(Coor[j](t), t)) then C[j] := i-1; V[j] := i;

print(Coor[j], C[j], " ", diff(Coor[j](t), t), V[j])

fi;

od;

od;

end:

>dx := diff(x(t), t, t) = -G*Mz*x(t)/(x(t)^2+y(t)^2)^(3/2);

>dy := diff(y(t), t, t) = -G*Mz*y(t)/(x(t)^2+y(t)^2)^(3/2);

>IniC := x(0) = 7*10^6, (D(x))(0) = 1, y(0) = 0, (D(y))(0) = 9*10^3;

>Ns := dsolve({dx, dy, IniC}, {x(t), y(t)}, numeric);

>dsnumsort(Ns(0), [x, y]):

>for i from 0 to 328 do T := 40*i; NsT := Ns(T);

>X[i], Y[i], Vx[i], Vy[i] := (eval([x(t), y(t), diff(x(t), t), diff(y(t), t)], NsT))[]

>od:

>i := 'i'; j := i+1;

>with(plots);

> plot([seq([X[i]-XS[j], Y[i]-YS[j]], i = 0 .. 328)], labels = ['dx', 'dy']);

>plot([seq([Vx[i]-VxS[j], Vy[i]-VyS[j]], i = 0 .. 328)], labels = ['dx', 'dy'])

> R := sqrt(XS[j]^2+YS[j]^2);

>sin := YS[j]/R; cos := X[j]/R;

> GCX := (X[i]-XS[j])*cos-(Y[i]-YS[i])*sin;

> GCY := (X[i]-XS[j])*sin+(Y[i]-YS[i])*cos;

> GCVx := (Vx[i]-VxS[j])*cos-(Vy[i]-VyS[i])*sin;

> GCVy := (Vx[i]-VxS[j])*sin+(Vy[i]-VyS[i])*cos;

> plot([seq([GCX, GCY], i = 0 .. 328)], labels = ['dx', 'dy']);

>plot([seq([GCVx, GCVy], i = 0 .. 327)], labels = ['dVx', 'dVy'])