restart;
with(LinearAlgebra): #Don't use the old package linalg
Q := x->0.5*erfc(x/sqrt(2));
rho:=0.25;sigma := 0.3;
m1 := sigma^2*Matrix(3, 3, [1, rho, rho^2, rho, 1, rho, rho^2, rho, 1]);
SNR := 15; mu := evalf(10^(.1*SNR));
#Changed a to a row vector:
a := <(ln(y/mu)+4*sigma^2)^2 | (ln(y1/mu)+4*sigma^2)^2 | (ln(y2/mu)+4*sigma^2)^2>;
b := MatrixInverse(m1);
#Changed c to a column vector
c := <(ln(y/mu)+4*sigma^2)^2, (ln(y1/mu)+4*sigma^2)^2,(ln(y2/mu)+4*sigma^2)^2>;
#Now a.b.c is a scalar, not a 1x1 matrix:
a.b.c;
L := exp(-(1/32)*a.b.c)/((32*Pi)^1.5*sqrt(Determinant(m1))*y*y1*y2);
#Using a more compact version of the integral:
BER:=evalf(Int(L*Q(sqrt((1/18)*y)+sqrt((1/18)*y1)+sqrt((1/18)*y2)),[y = 0 .. infinity, y1 = 0 .. infinity, y2 = 0 .. infinity]));
#The result comes pretty quickly (in Maple 16):
BER := 0.9862125481e-3

restart;
f:=sqrt(x^2+1)*ln(x^2+sqrt(x^2+1)); #Wrong integrand: See my other answer
A:=Int(f,x);
IntegrationTools:-Change(A,x=sinh(t));
simplify(%) assuming t>0;
Ft:=simplify(value(%)) assuming t::real;
Fx:=eval(Ft,t=arcsinh(x));

##Note: See comments below

F := sqrt(1+diff(y(x),x)^2);
Physics:-diff(F,diff(y(x),x));

The answer is not what you said, however.
Your answer is what you get by differentiating w.r.t. x as Carl notes.

Multiplying by cos(xx) to get rid of the singularities of tan helps.

I don't see eq1 and eq2. But supposing eq1 is the first argument to solve:
eq1:=(alpha+(l+alpha)*u+alpha*k*u^2)*a=
u*(alpha+(l+alpha)*u+alpha*k*u^2)*(1+l*alpha*b/((alpha+(l+alpha)*u+alpha*k*u^2)));
#Then you could use fsolve:
params:={l=1,alpha=1,b=100,k=50};
points:=seq([a,fsolve(eval(eq1,params),u)],a=0..20);
#Let eq2 be:
eq2:=f(a,u)=g(a,u); #Since I don't know what it is.
E:=map(x->[a=x[1],u=x[2]],[points]);
subs~(E,eq2);
#Or if you prefer:
eval~(eq2,E);

Of course your eq1 can be written simpler.

Try inserting two print statements in your procedure 'solution' just after the statement X[1,n]:=X[1,n]: (which statement doesn't do anything: Did you mean to unassign X[1,n]?).

print(f[1]);#################   
X[2,n] := diff(X[1,n], theta) - f[1]: ####RECURSIVE
print([2,n],X[2,n]);###############

After that run your program. You will see that f[1] includes the name X[2,1], and in the call solution(1) you have n=1, so you are making a recursive assignment where indicated above.

You may be interested in trying the parameters option:

solpar:=dsolve(eq1 union ICs,numeric,parameters=[A,N_h,b,ga,m,bi_h,bi_v,mu_h,mu_v]);
#Setting specific parameters:
solpar(parameters=[A = 5000., N_h = 10000., b = 1., ga = .167, m = 0., bi_h = .4, bi_v = .4, mu_h = 0.457e-4, mu_v = .25]);
#Plotting:
plots:-odeplot(solpar,[[t,I_h(t)],[t,I_v(t)],[t,R_h(t)],[t,S_h(t)],[t,S_v(t)]],0..50);
#Setting different parameters (only m is changed in this example):
solpar(parameters=[A = 5000., N_h = 10000., b = 1., ga = .167, m = 10000, bi_h = .4, bi_v = .4, mu_h = 0.457e-4, mu_v = .25]);
plots:-odeplot(solpar,[[t,I_h(t)],[t,I_v(t)],[t,R_h(t)],[t,S_h(t)],[t,S_v(t)]],0..50);

Let B be a matrix, not a list.
restart;
B := Matrix(1 .. 10, 1 .. 3);
for j to 3 do
  for i to 10 do
    B[i, j] := [RootFinding:-Analytic(eval(x*tan(x)-I*i+10), x = (j-1)*Pi+0*I .. j*Pi+10*I)]
  end do
end do:
#########
#Then to get values for i=1..10 and j=2 do
B[1..10,2];

The following works for your example:

f:=(x,y)->[x,`if`(y>0,4*y,y)];

or this weirder version:

f:=(x,y)->[`if`(x<1,x^2,x*y),`if`(y>0,4*y,y)];

Shouldn't Matlab:-FromMatlab be able to do it?

From the way you are using Psi I see that Psi is just a constant as are mu and c.
If we let u be the quantity obtained by evaluating (D[1, 1](f))(x1, x2, x3, x4) at the values in your worksheet, then u is rather large and unwieldy, even if you try evaluating it for Psi=1.
The question is: what do you hope to achieve? A smaller sized expression?
Is a series expansion in Psi of interest (about Psi=0) then
series(u,Psi=0,3);
If you want to expand all positive powers then that can be done
w:=evalindets(u,`^`(anything,posint),expand):

If we may assume that c is the speed of light then conceivably
asympt(u,c,3);

might be interesting.

From the help page for implicitplot3d:

"It uses an numerical algorithm based on triangulation into tetrahedrons combined with simple interpolation."

I guess the plotting interface is doing the interpolation and the resulting points are not available. That may be the reason that isosurface structures are not handled by plottools:-getdata.

You have assigned to tau[0], tau[1], and tau[2]. Thus implicitly tau has become a table. But you are also using tau as a name. This is known to cause serious problems. One solution is to change tau[0], tau[1], and tau[2] to tau0,tau1, tau2.
The same comment applies to your use of k.
Also: Change sum to add as you have concrete summations.
linalg is deprecated. LinearAlgebra should be used instead. But I don't see the need for linear algebra.

Finally, in A your F(tau) should be just F. The same goes for T,P, and Q.

restart;
S := Matrix( [ [a_1, a_2, a_3, a_4], [b_1, b_2, b_3, b_4], [c_1, c_2, c_3, c_4], [d_1, d_2, d_3, d_4] ] );
eq1:=a_1*d_1 = b_1*c_1, a_2*d_2 = b_2*c_2, a_3*d_3 = b_3*c_3, a_4*d_4 = b_4*c_4;
R := Matrix( [ [s_1, t_1, r_1, l_1], [s_2, t_2, r_2, l_2], [s_3, t_3, r_3, l_3], [s_4, t_4, r_4, l_4] ] );
eq2:=s_1*l_1 = t_1*r_1,s_2*l_2 = t_2*r_2,s_3*l_3 = t_3*r_3,s_4*l_4 = t_4*r_4;
T := Matrix( [ [1, 0, 0, 1], [0, 1, 0, 0], [0, 1, 0, 0], [1, 0, 0, 1] ] );
S.R=~T;
convert(%,set);
res:=solve(% union {eq1,eq2});
nops([res]);
res[1]; # First of the ninety

DE := diff(x(t), t) = x(t);                   
DF := diff(y(t), t) = k*y(t); #You missed k

xres:=dsolve(DE);
yres:=subs(_C1=_C2,dsolve(DF));
eq:=abs(y(t))=C*abs(x(t))^k;
eval(eq,{xres,yres}) assuming real;
combine(expand(%)) assuming real;

So indeed there is a constant C s.t. eq is satisfied (if that was what you wanted to show?).