By using solve with symbolic a and b, you will get 4  "solutions". However, they are not solutions all of them, because solve squares both sides of the equation to get rid of the square root. By doing that it includes the solutions of the equation having opposite sign on one side of the equation.
When a and b are of type numeric, solve is more careful.

restart;
eq:=3*a*(1+z)^2 = 2*(a*(1+z)^3+b)^(1/2);
solve(eq,z);
#Notice the 4th degree polynomial inside RootOf
#Now we do manually as said above:
map(x->x^2,eq);
#We have our 4th degree polynomial (with z instead of _Z):
collect((lhs-rhs)(%),z);
#Now let a and b be concrete from the start:
eq1:=eval(eq,{a=1,b=2});
#solve finds 2 solutions this time (now notice the indices singling out roots of the polynomial):
solve(eq1,z);
convert([%],radical); #If you want to see.
evalf(%);

So one has to be careful.

If your system is linear as it is in your very simple example, then this kind of thing can be handled by gaussian elimination.

L:=[x+y=0,2*x+2*y=0];
with(LinearAlgebra):
T:=GenerateMatrix(L,[x,y],augmented);
GaussianElimination(T);
#A less trivial case (the 3rd equation is obtained by adding the two first):
L2:=[x+y=1,x+2*y=3,2*x+3*y=4];
T2:=GenerateMatrix(L2,[x,y],augmented);
GaussianElimination(T2);

First of all notice that this is only (?) done when solving odes numerically as will be seen below.

I think the basic reason for turning a higher order ode into a system is that of efficiency. One procedure handles all orders. In theoretical tratments this is also very often done.

To see the actual statement you can look at the help page
?dsolve,numeric,ivp
Now this page only speaks of initial value problems, but actually the same is done for boundary value problems as will be seen below.

restart;
ode2:=diff(x(t),t,t)+4*x(t)=0; #Example
ode3:=diff(x(t),t,t,t)+8*x(t)=0; #Example
infolevel[dsolve]:=2:
#Exact (symbolic) solutions:
dsolve(ode2);
dsolve(ode3);
#Notice the specific mention of 2nd and 3rd order. For higher orders this is not done: Try 4th order.
#Now numeric (ivp):
dsolve({ode2,x(0)=0,D(x)(0)=1},numeric);
dsolve({ode3,x(0)=0,D(x)(0)=1,(D@D)(x)(0)=0},numeric);
#Now numeric, but bvp:
debug(`dsolve/numeric/bvp`);
dsolve({ode2,x(0)=0,x(1)=2},numeric);
#In the midst of the lengthy printout you will see the first order system written out.
#It appears in fproc, e.g.
#fproc := proc (M, t, Y, L, F) F[1] := Y[2]; F[2] := (-1)*4.*Y[1]; 0. end proc
#You may want to finish with
undebug(`dsolve/numeric/bvp`);

Suppose you have found an acceptable initial vector X0 (see my other answer).
Then a numerical solution could be found like this:
X:= < seq(x[i](t),i=1..7 )>;
#The ode system
sys:=convert(diff~(X,t)=~A.X,set);
#The initial conditions:
ics:=convert(eval(X,t=0)=~X0,set);
#The numerical solution:
res:=dsolve(sys union ics,numeric);
plots:-odeplot(res,[seq([t,x[i](t)],i=1..7)],0..0.1);

You may want to try MatrixExponential because it is so seemingly simple:

restart;
A:=ImportMatrix("G:/MapleDiverse/MaplePrimes/A-1.txt"):
C:=ImportMatrix("G:/MapleDiverse/MaplePrimes/C.txt"):
Y0:=ImportMatrix("G:/MapleDiverse/MaplePrimes/Y0.txt"):
Digits:=15:
E:=LinearAlgebra:-MatrixExponential(A,t):
simplify(fnormal~(eval(E,t=0),8)); #Just a check: Should be the identity matrix
#Considering you previously posted question I assume you may accept the following as X(0) (?):
X0:=LinearAlgebra:-LeastSquares(C,convert(Y0,Vector));
#The solution:
X:=E.X0:
plot([seq(Re(X[i]),i=1..5)],t=0..0.1); #Just to see something.

No solution.

A:=ImportMatrix("G:/MapleDiverse/MaplePrimes/A.txt");
Y:=ImportMatrix("G:/MapleDiverse/MaplePrimes/Y.txt");
LinearAlgebra:-LinearSolve(A,Y);
interface(rtablesize=infinity);
LinearAlgebra:-GaussianElimination(<A|Y>);

I had a brief look at it. Make sure that floats are used all over. An evalf here and there will help.

It seems that you are solving a second order ode by a difference scheme. What is the point? Why not solve the ode numerically in Maple using dsolve/numeric? Unless the point is to see how well the difference scheme works for different stepsizes?

Torq:= t-> 1:
m := 1;
L := 1;
g := 9.81;
ode:=m*L^2*diff(x(t),t,t)=Torq(t)-m*g*L*sin(x(t));
res:=dsolve({ode,x(0)=0,D(x)(0)=0},numeric,range=0..12);
plots:-odeplot(res,[[t,x(t)],[t,diff(x(t),t)],[t,rhs(ode)]],0..12,refine=1);

This seems to be related to

http://www.mapleprimes.com/questions/142322-Why-Has-The-Plot-Command-Deteriorated#comment142326

And indeed, it is the same problem:

restart;
eqn3:=diff(y(x),x)=x+sin(y(x));
IC3:=[y(0)=2];
DEtools[DEplot](eqn3,[y(x)],x=-5..5,IC3,y=-5..5,arrows=none,linecolor=BLACK);
plottools:-getdata(%);
M:=%[-1];
has(M,HFloat(undefined));
interface(rtablesize=infinity);
M;
M1:=M[9..-12,..];
plot(M1);

I hope this won't have to wait for Maple 17!

Added: The following simple remedy works (at least in the present case):

restart;
eqn3:=diff(y(x),x)=x+sin(y(x));
IC3:=[y(0)=2,y(0)=3];
DEtools[DEplot](eqn3,y(x),x=-5..5,IC3,y=-5..5,linecolor=black); p:=%:
eval(p,[undefined,undefined]=NULL);

#A more general version also works for the case in the link:

evalindets(p,Or(Matrix,listlist),m->remove(hastype,convert(m,listlist),undefined));

@Axel Vogt Here is a version inspired by your worksheet. It introduces a parameter 'a' replacing the number 3. Differentiates w.r.t. a and then integrates w.r.t. a from 1 to 3:

restart;
Int((ln(sin(x)^2+1)-ln(cos(x)^2+1))/(-cos(2*x)), x = 0 .. Pi/2);
with(IntegrationTools):
#Integrating only to Pi/4:
Int((ln(sin(x)^2+1)-ln(cos(x)^2+1))/(-cos(2*x)), x = 0 .. Pi/4);
combine(%);
#Using symmetry, the original integral to Pi/2 is
2*Change(%,x=arccos(t)/2);
combine(%);
Ba:=eval(%,3=a);
diff(Ba,a);
dBa:=value(%) assuming a>1;
#We find
eval(Ba,a=1);
solve((t+1)/(1-t)=u,t);
Change(%%,t=%);
value(%);
int(dBa,a=1..3)+%;
# Output:   arctan((1/4)*sqrt(2))*Pi

I think you must integrate numerically. But N is huge where H is not and vice versa. So some scaling may help.

restart;
eta:=10^(-2.623+(3388.8/(T-491)));
DG:=51749.47800-33.65211624*T-0.4480891609E-2*T^2;
N:=((T/eta)*(exp(108.3606-(3.02397E13/(T*(DG^2))))));
plot(N,T=500..1307);
#I have removed abs from the expression B since it is unnecessary in the range considered:
solve(51749.478-33.6521*A-0.0044809*A^2,A);
B:=(1.5438E-6*A*(1307-A))/(10^(-2.623+3388.8/(A-491)))*(1-exp(-(51749.478-33.6521*A-0.0044809*A^2)/(8.31*A)));
G:=Int(B,A=500..T);
H:=G^3;
#Now scale the two factors:
N1:=1e-15*N;
H1:=simplify(1e15*H);
J1:=N1*H1;
q:=evalf[15](Int(J1,T=500..1307));
q:=evalf[20](Int(J1,T=500..1307)); #Takes much longer because software floats must be used
Rcc:=evalf[20]((((4*Pi/3)*q)/1E-6)^(1/4));

I tried the following in Maple 16 (no problem in Maple 15):

restart;
u:=sqrt(75-3*x^2);
plot(u); #No graph shown
dt:=plottools:-getdata(%);
M:=dt[-1];
M[1..10,1..2]; #OK data
M[-10..-1,1..2]; #Last data point undefined
M1:=M[1..-2,1..2]; #Last data point removed
plot(M1); #OK
plot(M); #No curve
plot(u,smartview=false); #Doesn't help
plot(u,style=point); #Works!
plot(u,x=-6..6); #No graph
plot(u,x=-5..5); #OK

I also tried:

plot(x^(-1/3)); #As expected and as in previous versions
dt:=plottools:-getdata(%);
M:=dt[-1];
has(M,HFloat(undefined)); #Answer 'false'
#The same question about the M above was 'true'.

We may assume that 3 ending parentheses are missing at the very end (?).

indets(R);
Rc:=collect(R,[x,y],distributed);
##[Length of output exceeds limit of 1000000]
##But who would want to see it?
degree(Rc,{x,y});
C:=coeffs(Rc,[x,y],'L');
L;
L[5],C[5];

Here are two ways to do the manipulation. The second in three versions.

restart;
A:=LinearAlgebra:-RandomMatrix(1021,2); #Example matrix (array)
f:=(x,y)->x+2*y; #Example function
E:=zip(f, A[..,1], A[..,2]);
E[1..5];
M:=convert(E,Matrix); #May be needed for the export.

#A more straightforward way in 3 versions:
#Here the first column is kept intact:
B1:=Matrix(1021,2,(i,j)->if j=2 then f(A[i,1],A[i,2]) else A[i,j] end if);
A[1..5,1..2];
B1[1..5,1..2];
#Here the two first columns are intact but a third added:
B2:=Matrix(1021,3,(i,j)->if j=3 then f(A[i,1],A[i,2]) else A[i,j] end if);
B2[1..5,1..3];
#And this one only contains the results:
B3:=Matrix(1021,1,(i,j)->f(A[i,1],A[i,2]));
B3[1..5,1];

What type do the entries have?

B disappeared from your equations. I shall assume that the X on the right hand side ought to be B.

I had no problems with the following:

with(LinearAlgebra):
m:=400:
A:=LinearAlgebra:-RandomMatrix(3*m+1,2*m,generator=-1.0..1.0,outputoptions=[datatype=float[8]]);
B:=LinearAlgebra:-RandomVector(3*m+1,generator=-1.0..1.0,outputoptions=[datatype=float[8]]);
res1:=LeastSquares(A,B);
M:=LinearAlgebra:-Transpose(A).A;
res2:=LinearSolve(M,Transpose(A).B);
Norm(res1-res2);
res3:=M^(-1).Transpose(A).B;
Norm(res1-res3);

Is this close to what you want:

F:=Matrix([[x^2,sin(x)],[ln(x),cos(x)]]);
diff~(F,x);
FF:=unapply(%,x);
FF(Pi);