with(linalg):
A := matrix([[0, 1], [0, 0]]):
x0 := vector([0, 0]):
C := (0.01*psi1^2 + psi2^2)^(1/2):
t0 := 0:
T := 3:
step := 0.05:

gradC := subs(psi1 = psi1(t), psi2 = psi2(t), vector([diff(C, psi1), diff(C, psi2)])):
sol := []:
for alpha from 0 to 2*3.1416 by step do
	dsys := {
		diff(x1(t), t) = A[1,1]*x1(t) + A[1,2]*x2(t) + gradC[1], x1(0) = x0[1],
		diff(x2(t), t) = A[2,1]*x1(t) + A[2,2]*x2(t) + gradC[2], x2(0) = x0[2],
		diff(psi1(t), t) = -A[1,1]*psi1(t) - A[2,1]*psi2(t), psi1(0) = cos(alpha),
		diff(psi2(t), t) = -A[1,2]*psi1(t) - A[2,2]*psi2(t), psi2(0) = sin(alpha)
		}:
	sol := [sol[], dsolve(dsys, numeric, output = listprocedure)]:
end do:

X := proc(alpha) global step, T:
	eval(x1(t), sol[round(alpha/step)])(T);
end proc:
Y := proc(alpha) global step, T:
	eval(x2(t), sol[round(alpha/step)])(T);
end proc:

plot([X, Y, 0..2*Pi], numpoints=360);