JohnS asked about substitution in an "array/matrix", not in a Matrix or a list. 
This complicates matters a bit due to last-name evaluation.
The array and matrix types are deprecated now, but they still exist, so I'll try
to answer the question that was asked.  Here are two ways:

> A:= array([x,y,z]);

> subs(x=a,y=b,z=c,eval(A));

  [a, b, c]

> map2(subs,[x=a,y=b,z=c],A);

  [a, b, c]

This might be more sophisticated than what you want, but there is a procedure parametrization in the algcurves package that will provide a rational parametrization of an algebraic curve when that exists.  For example:

> eq_C := y=x^2+y^2;
  algcurves[parametrization](lhs(eq_C)-rhs(eq_C), x, y, t);

[(-1+t^2)/(2+2*t^2), (1-2*t+t^2)/(2+2*t^2)]

The <maple> tag is acting up again: that should be

[(-1+t^2)/(2+2*t^2), (1-2*t+t^2)/(2+2*t^2)]

In other words, the curve can be parametrized as x = (-1+t^2)/(2+2*t^2),
y = (1-2*t+t^2)/(2+2*t^2) (except for the point [1/2, 1/2] which is the limit as t -> infinity).

Suppose the file is myfile.csv in directory "c:\temp", consisting of two columns of data.
 

> M := ImportMatrix("c:/temp/myfile.csv", source=csv);
  P := CurveFitting[LeastSquares](M, t, curve=a0 + a1*t + a2*t^2);

Then to visualize the points and the curve:
 

> tmin, tmax := (min,max)(M[1..-1,1]);
   with(plots):
   display([pointplot(M), plot(P, t=tmin .. tmax)]);

plot3d(..., shading=zhue);

 See ?plot3d,options for other options

Tools, Options, Interface

1) Change "Default format for new worksheets" from Worksheet to Document

2) Change "Open worksheets in:" from "New Tab" to "New Window"

Then click "Apply to Session" if you only want it for this session, or "Apply Globally" if you want the change to be permanent

3) I don't think this is possible, short of saving the tab to a file and opening the file in a new window.

It looks to me like the DAE solver in dsolve can only do initial value problems, not boundary value problems.

As for Question 2, I think a similar hint to the one I gave for Question 1 should help.

k1(t) = -1, so k2(t) is nonzero except maybe at one point, so the same is true of k3(t),
so z(t) = 0 ...

 By the way, please don't ask the same question multiple times.

> evalf(BesselJZeros(0, 1..10));

2.404825558, 5.520078110, 8.653727913, 11.79153444, 14.93091771, 18.07106397, 21.21163663, 24.35247153, 27.49347913, 30.63460647

In the first example, from equation (6) either g(t) or k2(t) is 0.  If k2(t) = 0, then (2) says k1(t) = 0.  But this can't happen on any interval, because of equation (1).  So it must be g(t) = 0.  I'll leave the rest to you.

Do you mean you want a text file containing the four Arrays P, P1, P2, P3?

It would be somewhat easier if these were 51 x 3 Matrices rather than Arrays of lists.

> for Q in [P,P1,P2,P3] do
       fprintf("c:/mypath/myfile.txt","%{c,}V\n", 
            Matrix(51,3,(i,j)->Q[i][j]))
  end do:
  close("c:/mypath/myfile.txt");

If f(x) is defined, say, for 0 <= x < p, then f(x - p*floor(x/p)) will be periodic with period p and equal to f(x) for 0<= x < p.

You wrote your function using braces {} and square brackets[], where I think you mean parentheses ().  {} is used in Maple to enclose a set, while [] is used to enclose a list or an index.  For grouping, or the arguments of a function, Maple just uses (). 

After changing all { and [ to ( and all } and ] to ) I get a rather complicated function.  I think it would be a miracle if that equation had an explicit solution that could be written in closed form.  Maple certainly can't find such a solution.

Hint: For real a, there are three real roots counting multiplicity, so "exactly two distinct real roots" would require one of the roots to have multiplicity 2.  The discriminant of the polynomial is 0 when that happens.
 

LinearAlgebra[MatrixInverse](..., method=pseudo);

The only options for the arrows in fieldplot3d are LINE, THIN, SLIM, THICK, and `3-D`.
An alternative might be to compute a list of individual vectors and plot them using either plottools[arrow] or plots[arrow].

Perhaps something like this?

> with(plots): with(plottools):
    circle1:= circle([0,0],3):
    circle2:= circle([2,0],2):
    F:= (x^2+y^2-3^2)*((x-2)^2+y^2-2^2):
    region := implicitplot(F<0,x=-3..4,y=-3..3,filledregions=true,gridrefine=3):
    display([circle1,circle2,region], scaling=constrained, axes=box);