Example:

n:=5:
p:=add(b[n-i]*z^i,i=0..n)/add(a[n-i]*z^i,i=0..n);
asympt(p,z);
#Same as:
taylor(subs(z=1/t,p),t=0);
subs(t=1/z,%);

You could use fdiff (but maybe not at x=0).
This is the heat equation taken from the help page for pdsolve, numeric:

PDE := diff(u(x,t),t)=1/10*diff(u(x,t),x,x);
IBC := {u(x,0)=1, u(0,t)=0, D[1](u)(1,t)=0};
pds := pdsolve(PDE,IBC,numeric);
res:=pds:-value(t=0.8,output=listprocedure); ## t = 0.8 chosen as an example
#Now the procedure computing u(x,0.8) values:
U1:=subs(res,u(x,t));
#The partial derivative of u(x,t) at x = 0.5 and t=0.8:
fdiff(U1,[1],[0.5]);

Your boundary conditions are given at 3 points: 0, v, and L. Was that intended?
You mention an attached file. I don't see it.
You have a second order ode. You can only have two boundary conditions. You have 4.
You need to solve numerically. All constants must be numerical.

Since you are not getting an analytical expression by the the method you seem to want, why not just do like this:

restart;
EQ := -(1/781250)*exp(-(1/1250)*X[1]^2)*X[2]*exp(-(1/1250)*X[2]^2)/((1+exp((1/67)*X[2]-(1/67)*X[1]))*Pi);
U:=unapply(Int(EQ,X[2]=-infinity..infinity,method=_d01amc),X[1]);
evalf(U(3));
plot(U,-200..200);

Set maxfun=0. Contrary to what one might expect, this means that maxfun=infinity.

You may try Minimize from the Optimization package, but be aware that in general only a local minimum is found.

A:=1/2*x^2*y^2 + y^2*z^2 + z^2 *x^2 + 96/(x + y + z + 1);
cstr:=x^2 + y^2 + z^2 = 5;
Optimization:-Minimize(A,{cstr},assume=nonnegative);
#Visual check:
Z:=solve(cstr,{z});
plot3d(subs(Z[1],A),y=0..sqrt(5-x^2),x=0..sqrt(5));

But compare with:

B:=subs(Z[1],A);
res:=solve({diff(B,x)=0,diff(B,y)=0},{x,y}):
evalf(%);
remove(has,[%],I);
eval~(B,%);

So yes, only a local minimum was found by Optimimization:-Minimize. You could also look at DirectSearch available from the Maple Application Center.

If you denote your equations eq1, eq2, ..., eq6, then you will see that

eq1-eq2;

returns
-41656.25000*a[2, 1](t)-(250000/3)*a[1, 1](t)+12.47395837*a[2, 2](t)+41.6145833*a[1, 2](t)+8337.797619*a[2, 3](t)^2+10.41666667*a[1, 3](t)^2-(25/3)*a[2, 3](t) = 0

Since this is not an identity, this means that solving for the second derivatives (as has to be done by dsolve/numeric via DEtools[convertsys]) is impossible.

There is a similar problem with eq4-eq5;

You could also try the little linear algebra involved in solving for the second derivatives:

sys:=[seq(cat(eq,i),i=1..6)]:
sys0:=subs([diff(a[1, 1](t), t, t)=d11,diff(a[1, 2](t), t, t)=d12,diff(a[1, 3](t), t, t)=d13, diff(a[2, 1](t), t, t)=d21, diff(a[2, 2](t), t, t)=d22,  diff(a[2, 3](t), t, t)=d23],sys);
A,RHS:=LinearAlgebra:-GenerateMatrix(sys0,[d11,d12,d13,d21,d22,d23]);
LinearAlgebra:-GaussianElimination(A);
The rank of A is 4, not 6.

Differentiate eq3. (I'm assuming that you meant sqrt(1-tau^2) ).

restart;
eq1:=diff(x(t),t)=-z(t);
eq2:=diff(y(t),t)=z(t)-x(t)/2;
eq3:=y(t)=3-2*Int((z(t)-z(tau))/sqrt(1-tau^2),tau=0..1);
eq3a:=IntegrationTools:-Expand(eq3); #Somewhat better than eq3 in this context
diff(eq3a,t);
value(%);
eq4:=eval(eq2,%);
#Now solve the planar linear system in x and z:
dsolve({eq1,eq4,x(0)=1,x(1)=0});
res:=combine(%);
#To find y:
Z:=unapply(subs(res,z(t)),t);
eval(eq3a,z=Z);
resy:=value(%);
#The integral is not computed, but it is independent of t and can be evaluated numerically.
ynum:=rhs(evalf(resy));
plot(ynum, t=0..1);

Try this:

a:=<1,2;3,4>;
b:=Matrix(2,2,[1,2,3,4]);
addressof(a);
addressof(b);
#You see that these two matrices are stored at different locations. Thus evalb will find them different.
#However:
LinearAlgebra:-Equal(a,b);

See

http://www.mapleprimes.com/questions/146153-Boundary-Conditions

Try removing 4 of the boundary conditions:
IBC := {ui(-1000, t) = 10^(-5), ui(1000, t) = 10^(-5), ui(x, 0) = 0, ur(-1000, t) = 10^(-5), ur(1000, t) = 10^(-5), ur(x, 0) = 1/cosh(x)};
sol := pdsolve({l1, l2}, IBC, funcs, numeric, time = t);
sol:-plot([ui,ur],t=10);
sol:-animate([ui,ur],t=500,frames=100);

Set printlevel to 2 (or above):

printlevel:=2:

or use an explicit print statement.

Consider a simple example, in 3 versions:

restart;
for i from 1 to 5 do
 p:=i^2;
 if p=9 then "HaHaHa" end if
end do;
printlevel:=2:
for i from 1 to 5 do
 p:=i^2;
 if p=9 then "HaHaHa" end if
end do;
#Or with an explicit print statement and with printlevel 1
printlevel:=1:
for i from 1 to 5 do
 p:=i^2;
 if p=9 then print("HaHaHa") end if
end do;

Here is an example, where I use the DirectSearch package. It can be found in the application center.

restart;
ode:=diff(x(t),t)=sin(x(t))+sin(t^2);
res:=dsolve({ode,x(0)=0},numeric,output=listprocedure);
plots:-odeplot(res,[t,x(t)],0..6);
X:=subs(res,x(t));
#Assuming that DirectSearch is in the following location "G:/MapleDiverse/DirectSearch2/EnglishVersion"
libname:="G:/MapleDiverse/DirectSearch2/EnglishVersion",libname;
with(DirectSearch);
MI:=GlobalSearch(X(t),[t>=0 and t<=6]);
MA:=GlobalSearch(X(t),[t>=0 and t<=6],maximize);

Obviously the first one is wrong.
Take a very special case where the integral can be computed:

f:=int(y*exp(-y+I*y*omega), y = 0 .. infinity);
#Resulting in    f:=-1/((2*I)*omega+omega^2-1)
#Then we get
fourier(f,omega,t);
#resulting in  2*exp(-t)*Pi*t*Heaviside(t)

The variable in in a loop (in your case 'tend') takes on a value.
When you do
for tend from 0 to 50 do ..... end do;
the variable tend takes on the values 0, 1, 2, ... consecutively.
(If you didn't have the while clause tend would end up with the value 51. With the while clause it ends up being 13.)
Since at each iteration tend is a number not a name, tend=0..50 makes no sense.
Your third attempt has the same problem.
In your second attempt 'tend' is never given a value, so whether tend <= 12 or not cannot be determined.