The initial conditions for the derivatives must be given as follows:

D(x)(0) = 15, D(y)(0)= 1, D(z)(0) = 0

When that change has been made you get a solution and a plot of all 3 together can be obtained e.g. this way

plots:-odeplot(sol1,[[t,x(t)],[t,y(t)],[t,z(t)]],0..10);

With k just a name p(k) evaluates to 0 in Maple. Thus so does  sum(p(k), k = 1 .. 1), since sum obeys normal evaluation rules.

Delaying evaluation by using unevaluation quotes solves the problem.

sum('p(k)', k = 1 .. 1);

add on the other hand has special evalution rules (like seq does) so here no quotes are necessary:

add(p(k), k=1..1);

In your first loop you define theta[k+1], but since that involves theta[k+1] itself (for i = 0) you need theta[k+1] to be known already. But it isn't. To begin with only theta[0] is known. So when k = 0 you are attempting to define theta[1] from theta[1], which leads to the infinite recursion error.

eq:=diff(u(r,t),t)=diff(u(r,t),r,r)+diff(u(r,t),r)/r+u(r,t)+f(r,t);

##
#Picking as an example f(r,t)=r^2*t
##

res:=pdsolve(eval(eq,f(r,t)=r^2*t),{u(r,0)=0,u(0,t)=t^2,u(2,t)=0},numeric);
res:-plot(t=1,title="The solution at time t = 1");
res:-animate(r=0..2,t=0..3,title="t=%f");
res:-animate(t=0..3,r=0..2,title="r=%f");

If f is defined as a function ( and if I understand your intention correctly) you could use D as in this example,

f:=sin:
plot(D(f)(2*Pi/z),z=1/(2*Pi)...60/(2*Pi)); # plots the derivative

plot(f(2*Pi/z),z=1/(2*Pi)...60/(2*Pi)); # plots the function itself

The plot looks like the plots produced by the command line version of Maple, cmaple.

I have seen this phenomenon occur in Maple 14 on a few rare occasions for no apparent reason. I can't reproduce it though, since I have no clue why it happens when it does. I tried your code and the result was just fine in both Maple 14 and 15.

Your code stripped of output:

with(plots);
homo_data:= [[0.07143, -5.27798], [0.03571, -5.09811],[0.02381, -5.04314], [0.01786,-4.95879]];
lumo_data:= [[0.07143, -2.0953], [0.03571, -2.52144],[0.02381, -2.62239], [0.01786,-2.79138]];
homo_eqn:= -12.94151*sqrt(1-0.85175*cos(Pi/(N+1)));
eval(homo_eqn, N=infinity);
lumo_eqn:=homo_eqn+15.90268*sqrt(1-0.98083*cos(Pi/(N+1)));
gap_data:=[[0.07143, 2.9816], [0.03571, 2.3132],[0.02381, 2.1291], [0.01786, 2.0525]];
display([pointplot(homo_data),pointplot(lumo_data),plot([[1/N,homo_eqn(N),N=14..10000],[1/N,lumo_eqn(N),N=14..10000]], labels = ["1/N","E in eV"])]);

Try those lines in a new worksheet and see what happens.

After option package;  insert the line

uses plots;

Then all procedures in the plots package can be used with their short names, e.g.  display instead of plots:-display.

(A comment: The redefinition of `plots/animate` is not necessary in Maple 15, but is in earlier versions).

(rhs@op)~([a]);

For concrete values of m and n it is no problem, so if that suffices you don't have to do much:

ode:=diff(f(x),x,x)+f(x)^n=0;
f:=x->sum(u[i](x)*p^i,i=0..m);
ode;
eval(ode,{m=3,n=2});
collect(%,p);

Does this happen after a restart?

And a comment: If you want to plot sin on the interval -10..10 the syntax is

plot(sin(x), x=-10..10);

You can animate animate itself.

A not very exciting example:

animate(animate,[plot,[b*x^a,x=0..1,caption=typeset("a = ",a,", b = ",b)],a=0..5,frames=10,paraminfo=false],b=1..2,paraminfo=false,frames=5);

You have too many initial and boundary conditions. You are allowed 4 (1 + 1 +2, adding the orders of the 3 odes).

Also 'a' needs a value. Initial or boundary conditions at x1 for the derivative should be given using D(F)(x1)= whatever.

Example using only initial conditions and using a = 1:

dsys:={DN,Du,DDF,F(x1)=f1,u(x1)=u1,N(x1)=n1,D(F)(x1)=f_1}; 

dsoln:=dsolve(dsys,numeric);

plots:-odeplot(dsoln,[x,F(x)],x=x1..x2);

Try

convert(arctanh(x),ln);

then you will get 1/2*ln(x+1)-1/2*ln(1-x). 

For -1 < x < 1 both logaritms are real. For x < -1 the first logarithm ln(x+1) is the logarithm of a negative number and for x > 1 the second logarithm ln(1-x) is.

Maple uses the principal value of the complex logarithm which is defined as

ln(z) = ln(abs(z)) + I*argument(z)

where argument(z) is the principal argument of z, i.e. the argument v satisfying -Pi < v <= Pi.

Thus ln(x+1) = ln( abs(x+1)) + I*argument(x+1) = ln( abs(x+1)) +I*Pi, when x < -1.

Similarly, for x > 1, you get

ln(1-x) = ln( abs(1-x)) +I*Pi.

No evaluation takes place when defining a procedure (function). Thus in your loop q is not evaluated to the 3 different values. If you use unapply(p*q*x, x)  then p*q*x is evaluated and then the function is returned.

phi:=Array(1..3):

for q to 3 do

     phi[q] := unapply(p*q*x,x);

end do;

If the singularity here means that the solution or one of its components tends to infinity as the independent variable tends to the singularity, then continuing the solution on the other side would not make any sense. How should it be continued?

An example:

ode:=diff(x(t),t)=1+x(t)^2;
res:=dsolve({ode,x(0)=0});

The solution is x(t) = tan(t). It is correct that tan(t) is also defined for t>Pi/2, but why choose tan(t) for t>Pi/2, why not chose any other solution x(t) = tan(t+C) ?