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) ?

If you give us the code for euclid we could try the loop ourselves.

Here are two solutions. Both use unapply in the definition of f in order that the variable a in point21 and point22 is seen as the same as the formal variable a.

The first version corrects an old problem in `plots/animate`. The problem is that the occurrences of subs in that procedure ought to be replaced by eval. I hope that problem has been removed in Maple 15. We shall see.

restart;
S := solve([y-a^2 = -(x-a)/(2*a), y = x^2+1], [x, y]);
point1 := [a, a^2];
point21:=eval([x,y],S[1]);
point22:=eval([x,y],S[2]);
f := unapply(piecewise(abs(a)>=1/2, point22, point21),a);
p1 := plot(x^2, x = -2 .. 2, color = green);
p2 := plot(x^2+1, x = -2 .. 2, color = blue, scaling = constrained);
###Correcting the above mentioned problem:
`plots/animate`:=subs(subs=((x,y)->eval(y,x)),eval(`plots/animate`)):
###
anm := plots:-animate(plot, [[point1, f(a)]], a = -2 .. 2);
plots:-display(p1, p2, anm);

##Second version, which also works with the uncorrected animate

#Create a procedure producing just a plot of the line:
pp:=a->plot([[a, a^2], f(a)]);
#Now animate the plot procedure pp:
anm:=plots:-animate(pp, [a], a = -2 .. 2);
plots:-display(p1, p2, anm);

There are two roots. For a between -1/2 and 1/2 you must use S[1].

The result of

5*3 mod 7;

is 1, meaning that the multiplicative inverse of 3 is 5, when you calculate mod 7.

The command

solve(Prob(w,T)=0,w,AllSolutions);

gives you a formula for all solutions. Variables of the form _Z1~ represent integers, variables of the form _B1~ take values 0 or 1. Since sin is squared all the roots have multiplicity 2.

If your Eqs are meant to be equations, then b is an unknown too.

Eqs := [x*b = y*b, b = 0];
                      
solution := solve(Eqs, [x,y,b]);

I haven't examined your worksheet except for syntax and the like. But certainly T1 needs to have a numeric value if the odes are solved numerically.

You can use the parameters option, if T1 cannot be determined from whatever else is in your problem.

try1 := dsolve(eqnew union {ics}, numeric, parameters = [T1]);
try1(parameters=[1.2345]);
odeplot(try1,[t,a(t)],0..0.2);

Although the parameter is called t, in fact the parameter is x. If you make both x and y appropriate functions of time t then you can obtain a reasonable fall and rise of the ball.

The last formal parameter in your procedure Lplot is dl. When you call the procedure with dl as a sequence (as you did: the global dl is a sequence), only the first element is passed to the procedure.

You should consider the formal dl a set and pass a  set, like this: 

Lplot := proc (t0, xt, yt, zt, para, R, dl) global L;
data(t0, xt, yt, zt, para);
L := dsolve(`union`(subs({op(para)}, dl), {op(initial)}), {x(t), y(t), z(t)}, type = numeric, output = listprocedure, range = R)
end proc;
Lplot(0, 1, 1, 1, para, 0 .. 30, {dl});

PT:=LinearAlgebra:-RandomMatrix(3,generator=-9.0..9.);

X0:=<1.,2,3>:
X1:=<4,5,6>:

to 5 while LinearAlgebra:-Norm(X1-X0,2) > 0.0000001 do
X0:=X1;
X1 := PT.X1;
end do;


Or to handle the actual problem:

X0:=<1.,2,3>:
X1:=PT.X0:
for i to 1000 while LinearAlgebra:-Norm(X1-X0,2) > 0.0000001 do
X0:=X1/LinearAlgebra:-Norm(X1,2);
X1 := PT.X0;
print(`/`~(X1,X0));
X1:=X1/LinearAlgebra:-Norm(X1,2)
end do:
i;

The line

eq := z-> H^2=(-1/(12*f1))*(T*Omega+f):

needs to be replaced by

eq :=unapply( H^2=(-1/(12*f1))*(T*Omega+f), z):

or else eq(1.2345) would contain the global z present in Omega. You defined Omega by

Omega := H0^2*Omega0*(1+z)^3/H^2: