@Carl Love

pde1:=diff(u(x,t),t)=u(x,t); #Formally a pde
pdsolve({pde1,D[2](u)(x,0)=f(x)});
pde2:=diff(u(x,t),t)+diff(u(x,t),x,x)=0;
pdsolve({pde2,D[1,1](u)(x,0)=sin(x)});
pdsolve({pde2,D[1,1](u)(Pi/2,t)=exp(t),D[1](u)(Pi/2,t)=0});

@aybud
Actually the following simpler version also works (now Maple 17):

restart;
mystery := proc (z) (1+I)*z end proc;
F1 := proc (N, t, Y, YP)
   YP[1] := conjugate(mystery(Y[1]));
end proc;
curve := dsolve(numeric,procedure = F1, number = 1, procvars = [z(t)], start = 0, initial = Array([-.2+I*.3]));
plots:-odeplot(curve,[Re(z(t)),Im(z(t))],-1..1);

########################
But the next version, where mystery is made a function of 2 real variables works also with DEplot:

restart;
mystery := proc (x,y) local z; z:=x+I*y; (1+I)*z end proc;
F1 := proc (N, t, Y, YP) local tmp;
   tmp:=mystery(Y[1],Y[2]);
   YP[1] := Re(tmp);
   YP[2] := -Im(tmp);
end proc;
curve := dsolve(numeric,procedure = F1, number = 2, procvars = [x(t), y(t)], start = 0, initial = Array([-.2,.3]));
plots:-odeplot(curve,[x(t),y(t)],-1..1);
curve(1.2345); #No 0.*I appearing
inits := [seq(seq([0, -.2+.4*k, .3*j], k = 0 .. 1), j = 1 .. 4)];
DEtools[DEplot](F1,[x(t),y(t)],t=-1..1,number=2,inits);

@aybud
Actually the following simpler version also works (now Maple 17):

restart;
mystery := proc (z) (1+I)*z end proc;
F1 := proc (N, t, Y, YP)
   YP[1] := conjugate(mystery(Y[1]));
end proc;
curve := dsolve(numeric,procedure = F1, number = 1, procvars = [z(t)], start = 0, initial = Array([-.2+I*.3]));
plots:-odeplot(curve,[Re(z(t)),Im(z(t))],-1..1);

########################
But the next version, where mystery is made a function of 2 real variables works also with DEplot:

restart;
mystery := proc (x,y) local z; z:=x+I*y; (1+I)*z end proc;
F1 := proc (N, t, Y, YP) local tmp;
   tmp:=mystery(Y[1],Y[2]);
   YP[1] := Re(tmp);
   YP[2] := -Im(tmp);
end proc;
curve := dsolve(numeric,procedure = F1, number = 2, procvars = [x(t), y(t)], start = 0, initial = Array([-.2,.3]));
plots:-odeplot(curve,[x(t),y(t)],-1..1);
curve(1.2345); #No 0.*I appearing
inits := [seq(seq([0, -.2+.4*k, .3*j], k = 0 .. 1), j = 1 .. 4)];
DEtools[DEplot](F1,[x(t),y(t)],t=-1..1,number=2,inits);

@Carl Love Thanks, but what made me think something went before was among other things:
for i from low by byy to high do
proc (X) options operator, arrow; X(t) end proc;
proc (Y) options operator, arrow; Y(t) end proc;
.....
What is the purpose of the last two lines? And where are VXX, VYY, Cordx, Cordy defined?

Apparently you have 18 equations. Could you give us an idea of what they are all about, preferably by giving us a shorter example, alternatively by uploading a worksheet with the full set of equations?
"output=listprocedure" is an option for dsolve/numeric, not for fsolve.

Now what went before causing you to write "so i got this"?

How do you make sense out of that question, or is it one? "So i got rhis ... "

@LEETZ You can use value:

uval:=subs(s:-value(output=listprocedure),u(r,t));
#The first command gives the same as s:-animate:
plots:-animate(plot,[uval(r,t),r=0..1],t=0..5);
#Without using cylindrical coordinates:
plots:-animate(plot3d,[uval(sqrt(x^2+y^2),t),x=-1..1,y=-sqrt(1-x^2)..sqrt(1-x^2)],t=0..5);
#Using cylindrical coordinates:
addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]);
plots:-animate(plot3d,[uval(r,t),r=0..1,theta=0..2*Pi,coords=z_cylindrical],t=0..5);

@LEETZ You can use value:

uval:=subs(s:-value(output=listprocedure),u(r,t));
#The first command gives the same as s:-animate:
plots:-animate(plot,[uval(r,t),r=0..1],t=0..5);
#Without using cylindrical coordinates:
plots:-animate(plot3d,[uval(sqrt(x^2+y^2),t),x=-1..1,y=-sqrt(1-x^2)..sqrt(1-x^2)],t=0..5);
#Using cylindrical coordinates:
addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]);
plots:-animate(plot3d,[uval(r,t),r=0..1,theta=0..2*Pi,coords=z_cylindrical],t=0..5);

@geri23 You may want to take a look at the help page:

?dsolve,events

Here is a simple version of a  thermostat model inspired by Joe Riel a few years ago:

dsys := {diff(T(t),t) = piecewise(b(t)=0, -k*(T(t)-Tu), -k*(T(t)-Tu)+F)
          , T(0) = T0
          , b(0) = 1
        };

L := dsolve(dsys, numeric
                , 'events' = [[T(t)=T1, b(t)=1],[T(t)=T2, b(t)=0]]
                , 'discrete_variables' = [ b(t) :: boolean ]
                , parameters=[Tu,T1,T2,T0,k,F]):
L(parameters=[10,18,22,7,1,13]);
plots:-odeplot(L, [[t,T(t)]], 0..10);

@geri23 If that happens then it likely illustrates that discrete methods may be unstable under certain conditions, in particular if the stepsize is too large.
Why do you insist on a fixed stepsize method. What is the point? Is it perhaps to exhibit the instability?

@geri23 Apparently your time unit is second (T is 24 hours). So you may want to solve the system using Euler's method with stepsize = 5*60 (5 minutes).

So you don't insist that the matrix should be a rotation matrix, but instead that the sum of the rows and the first two columns all be 1 and the determinant be 1?

with(LinearAlgebra):
R:=Matrix([[R11,R12,fx],[R21,R22,fy],[R31,R32,fz]]);
Determinant(R);
Transpose(R).R=IdentityMatrix(3); #You don't require this, or do you?

Why the second root? The ordering of the output may not be what you want:

for alpha from .4 by .1 to 1 do
   argument~([solve(x^4-alpha, x)]);
end do;

Number 2 root has argument Pi/2 in all but the last case, where it is Pi.

You could sort the output e.g. in one of these ways:
for alpha from .4 by .1 to 1 do
   sort([solve(x^4-alpha, x)],(x,y)->evalb(Re(x)<Re(y) or Re(x)=Re(y) and Im(x)<Im(y)));
end do;
#or
for alpha from .4 by .1 to 1 do
   sort([solve(x^4-alpha, x)],(x,y)->evalb(argument(x)<argument(y)));
end do;

The advantage of Markiyan's version is that the output is the sequence sought. A print statement just prints to the screen. You have to copy the result from the screen afterwards.
Using a loop the same could be done with a small change:

S:=NULL:
for alpha from .4 by .1 to 5 do
   S := S,solve(x^4-alpha, x)[2];
end do:

S;