You may try this (done in Maple 12, as it so happens).

sol := dsolve(Systeme, type = numeric,range=0..20,output=listprocedure);

 odeplot(sol, [T1(t), T2(t)], 0 .. 20,refine=1);

TT1,TT2:=op(subs(sol,[T1(t),T2(t)]));
p:=t->pointplot([cos(TT1(t)), sin(TT1(t))],color=red,symbol = solidcircle, symbolsize = 25):
animate(p,[t],t=0..8,frames=300);

# Added later: In Maple 15 (as I suspected) you can just do

animate(pointplot,[[cos(TT1(t)), sin(TT1(t))],color=red,symbol = solidcircle, symbolsize = 25],t=0..8,frames=300);

#The need for defining the procedure 'p' above is due to a weakness in 'animate' present in Maple versions earlier than 15, but corrected in Maple 15.

Replace Epr11 := evalc(conjugate(psi1*psi1));

by

Epr11 := evalc(  conjugate(psi1(x))   *  psi1(x)  );

and similarly for the other quantities.

You may try to isolate the time derivatives:

SYS:=solve(sys,{diff(pA(z,t),t),diff(pB(z,t),t),diff(pC(z,t),t)}):
#However, the system is then rather huge:
length(SYS);
#You don't get the error message you got before, but the computation takes longer than my patience allows for:

pdsolve(SYS,ibc,[pA,pB,pC],numeric,time=t,range=0..1);

In defining 'ode' in your corrected version there is no (implied) multiplication sign after R in ode[2].

In defining B there is a factor 2 missing.

I would never use 2D math input, especially for serious work like this.

From what you are saying it seems that you are assigning to y(t). If that is the case, don't, but assign the right hand side of the output from dsolve to a name like 'yt' or 'Y'.

RealDomain has its own 'limit' (try replacing ':' by ';' in with(RealDomain) ) and that apparently has a bug.

However, the equivalent version

limit((1-x)^(-1/x),x=0,left);

gives the correct result.

By debugging the procedure RealDomain:-limit with your original version limit((1-1/x)^(-x),x=-infinity), we see that the local variable 'ball' gets assigned the range

RealRange(Open(Float(-infinity)), Open(Float(-infinity)))

and that may be the problem as it results in the local variable 's' being assigned the value 'false' because the statement

is((1-1/x)^(-x),real) assuming x::RealRange(Open(Float(-infinity)), Open(Float(-infinity)));

results in 'false'.

[ Notice that the statement
RealRange(Open(Float(-infinity)), Open(Float(-infinity)));
evaluates to RealRange(-infinity, Open(Float(-infinity)))
thus the range is of the form [a,a), which is meaningless if a is finite. But a = -infinity, so? ]

Finally that s = false makes the procedure return 'undefined'.

How about something like

dsol1 :=dsolve({dsys1,ini1},numeric,method=lsode,var, abserr=1e-9, relerr=1e-8,output=Array([seq(k*th,k=0..N1)]));

And the matrix of values:

dsol1[2,1];

To see whether there is any justification for this method you could look at a much simpler problem.

Here I have taken a very simple problem, where the exact solution is easily found, so that the actual error can be found. If there is no connection in this simple case, there is no reason to expect a connection in a more complicated case.

restart;
Digits:=15:
sys:={diff(y(t),t) = y(t), y(0) = 1};
##The exact solution is just the exponential function exp:
dsolve(sys);
#Solving numerically:
res:=dsolve(sys,numeric, method =rkf45,abserr=10^(-10),relerr=10^(-10),output=listprocedure);
Y:=subs(res,y(t));
#Shortening your definition of S1 some (not relevant):
N := 1000:
S1 := seq([0.005 * k, Y(0.005 * k)-Y(0.005 * (k-1))], k = 2..N):
#Plotting S1 and the actual error in the same system:
plot([[S1],exp(t)-Y(t)],t=0..5,thickness=3,legend=["S1","Actual error"]);
#Plotting the actual error only:
plot(exp(t)-Y(t),t=0..5,legend="Actual error");

#You could try other simple problems. Conclusion: The answer is No!

In fact all we need to realize is that S1 would not be a sequence of zeros even if the exact solution was used for its computation, i.e. even if a perfect method was used!

I haven't checked whether your assertion is right or wrong, but you could use the StringTools package like this:

showstat(IntegralTransforms:-Tools:-SimplifyPower);
STR:=convert(eval(IntegralTransforms:-Tools:-SimplifyPower),string):
StringTools:-Search("j := evalf(j[-1][2]-1);",STR);
StringTools:-Search("k := evalf(j[1][1]+1)",STR);
StringTools:-Search("j := evalf(j[-1][2]-1); k := evalf(j[1][1]+1)",STR);

STRnew:=StringTools:-Substitute(STR,"j := evalf(j[-1][2]-1); k := evalf(j[1][1]+1)","j := evalf(t1[-1][2]-1); k := evalf(t1[1][1]+1)"):
P:=parse(STRnew):
showstat(P);

Then overwrite (if you really want to) IntegralTransforms:-Tools:-SimplifyPower. (I didn't try).

There is a syntax error. 'sys' is a set and 'values' is a sequence. Thus you need

dsolve(sys union {values}, .....);

After that you find that the implicit option is not available.

You could try a numerical solution.

indets(sys,name);
res:=dsolve(eval(sys,{seq(l[k]=1,k=1..3),seq(m[k]=1,k=1..3),g=10}) union {values},numeric,output=listprocedure);
plots:-odeplot(res,[seq([t,theta[k](t)],k=1..3)],0..0.338);

The warnings tell you the problem. The solutions (edited: their derivatives, not their values) with initial values ivp3 become infinite in finite time (in fact in less than 7 time units). Adding a range to y (and a small stepsize) may help as seen below, but notice that warnings are not errors, just warnings.

g:=y->(800-y)*y/(100+y)-450;
solve(g(y)singular(g(y),y);
sys1:=diff(y(t),t)=g(y(t));
with(DETools):

ivs3:=[y(0)=100,y(0)=250,y(0)=800]:

DEplot(sys1,y(t),t=0..10,ivs3,y=0..800,stepsize=.01);

#An alternative approach is to consider instead dt/dy = 1/(right hand side of sys1):
sys2:=diff(t(y),y)=1/g(y);
ivs2:=[t(100)=0,t(250)=0,t(800)=0];
DEplot(sys2,t(y),y=-100..800,ivs2,t=0..10,arrows=line);
p1:=%:
#Now transform the plot:
T:=plottools:-transform((y,t)->[t,y]):
T(p1);

The culprit seems to be 'radialstart', which apparently has problems with discrete data (not only list of lists of reals).

restart;
with(plots):
g:=t->3+4*cos(t);
polarplot(g(theta),theta=0..2*Pi, radialstart=-2);
N:=5:
t:=[seq(2*Pi*(i-1)/(N-1),i=1..N)]:
L:=map(s->[g(s),s],t):
polarplot(L);
plottools:-getdata(%);
polarplot(L,radialstart=-2);
plottools:-getdata(%);

#A workaround could be:

L2:=map(s->[g(s)+2,s],t):
polarplot(L2,tickmarks=[[seq(k=k-2,k=1..8)],default]);

You could use animate, e.g. like in these examples.

animate(implicitplot,[(x^2+y^2-1)^3-x^2*y^3 = 0, x = -1.2 .. 1.2, y = -1.2 .. d, grid=[50,50],gridrefine=2],d=-1.1..1.3,paraminfo=false);
animate(implicitplot,[(x^2+y^2-1)^3-d*x^2*y^3 = 0, x = -1.2 .. 1.2, y = -1.2 .. 1.3, grid=[50,50],gridrefine=2],d=-1..1,paraminfo=false);
animate(implicitplot,[(x^2+y^2-d)^3-x^2*y^3 = 0, x = -1.2 .. 1.2, y = -1.2 .. 1.3, grid=[50,50],gridrefine=2],d=0..1,paraminfo=false);
animate(implicitplot,[(x^2+y^2-d)^3-x^2*y^3 = 0, x = -1.2 .. 1.2, y = -1.2 .. 1.3, grid=[50,50],gridrefine=2,filledregions,coloring=[red,white],color=white,axes=none],d=0..1,paraminfo=false);

A text search in Maple for 'mutable' results in (among other things) 'Programming Guide, Chapter 4'.

Here is a numerical solution. The "matrix singular" situation is avoided basically by replacing l^4 by the name l4, but also by introducing l4 as a constant function:

eqn2 := diff(Y(x), x$4) = l4(x)*Y(x),diff(l4(x),x)=0;
bc := Y(0) = 0, D(Y)(0) = 1, D(D(Y))(0) = 0, Y(.5) = 0, D(Y)(.5) = 0;
sol2 := dsolve({bc, eqn2},{Y(x),l4(x)},numeric,output=listprocedure);
plots:-odeplot(sol2,[x,Y(x)]);
plots:-odeplot(sol2,[x,l4(x)^(1/4)]);
L2:=subs(sol2,l4(x)^(1/4));
L2(0);
L2(.5);