It appears that you are using 2D math input and have chosen to define f as a function.

You might just as well have defined your expression as just that: an expression in x, like this:
f1:=2*sin(0.05*Pi*x-Pi*0.5)+2;

If you don't give a range, plot will choose one for you, sometimes -10..10 under some circumstances -2*Pi..2*Pi.
Try e.g.
plot(cos);
In your case do
plot(f, 0..40);
Alternatively:
plot(f(x), x = 0..40);
###
In the expression case (using f1) just do
plot(f1, x = 0..40);
###
The rewriting you are talking about must refer to the default treatment of Pi in products with a float, but that is irrelevant for your present plotting problem.
Try
Pi*0.5;
If you get 1.570796327 you have the default setting of floatPi=true in kernelopts.
That can be set to false, as I myself have done.
kernelopts(floatPi=false);
You can put that statement into your maple.ini file, if you wish.

You say that exp(alpha1*(Dp*alpha1*t-lh+x))  is infinite in one interval. That could only be the case if alpha1*(Dp*alpha1*t-lh+x) is infinite in one interval. But that quantity is finite for all t >=0 (and indeed all t, but here only positive t are interesting).

(You have a thetac and a Dc in your assumptions, but they are not in the expression.)

Basically what you have is:
 

eq1:=exp(-a*sqrt(s))/(sqrt(s)+b);
inttrans[invlaplace](eq1,s,t) assuming a>0,b>0;

which returns
exp(-(1/4)*a^2/t)/sqrt(Pi*t)-b*exp(b*(b*t+a))*erfc((1/2)*(2*b*t+a)/sqrt(t))

What you have in your definition of A and B are sequences. A:=(-4,5,-2) defines a sequence. You don't even need the parenthesis.
What you need to do is to define vectors e.g. by using <... > as in:
 

A,B:=<-4,5,-2>,<8,2,7>;
A-B;

There are several syntactical errors:
1. While odes is a sequence, ics is a set, so you must change the dsolve command, e.g. to dsolve({odes} union ics, .... ).
2. Did you intend to set Digits:=50? in that case change digits:=50 to that.
3. What was the intended meaning of the 4th ode:
 

ode4:=diff(r^4*exp(-v(r))*sqrt(1-2*G*M(r)/(r*c^2))*(diff(1-omega(r)/Omega, r)), r)+4*r^3*(1-omega(r)/Omega, r)*(diff(exp(-v(r))*sqrt(1-2*G*M(r)/(r*c^2)), r)) = 0;
## This part in particular:
op(5,lhs(ode4));
nops(%); # A product of 4 factors

The product of the first 3 factors is 4*r^3*(1-omega(r)/Omega, r). The last factor doesn't make any sense in this context.

You simply used = where := was called for in two places.
After correcting that, however, you get:

ExponentialFit(THAXexa, THAXexb, v);
Error, (in Statistics:-ExponentialFit) dependent values must evaluate to positive numbers

But then you could do
-ExponentialFit(THAXexa, -THAXexb, v);

####
Notice that THAXexa := THAX[[1 .. 10], [1]] defines a 10x1 matrix. If you want a vector do
THAXexa := THAX[1 .. 10, 1];
or simply
THAXexa := THAX[ .. , 1];
####
I forgot that you also wanted the datatype to be float. You could do
THAXexb := Vector(THAX[..,2],datatype=float);

You could do as follows:
 

sol:=dsolve(diff(y(x),x)= x/(sqrt(x^2-16))*1/(2*y(x)));
sol1:=subs(x^2=u+16,[sol]);
simplify(%) assuming u>0;
subs(u=x^2-16,%);

I managed to get your code from the link.
Here is a way of plotting. I'm using numerical solution:

restart;
with(plots):
##
beta[v] := 0.5e-2; beta[h] := 0.5e-2; alpha := .40; mu[h] := 1/(365*60); mu[v] := 1/14; omega := .60; Lambda[h] := 100; Lambda[v] := 1000; gamma1 := .9999167;
##
sist := (D(s[h]))(t) = Lambda[h]+omega*r[h](t)-(beta[h]*i[v](t)+mu[h])*s[h](t), (D(i[h]))(t) = beta[h]*s[h](t)*i[v](t)-(gamma1+mu[h])*i[h](t), (D(r[h]))(t) = gamma1*i[h](t)-(omega+alpha+mu[h])*r[h](t), (D(m[h]))(t) = alpha*r[h](t)-mu[h]*m[h](t), (D(s[v]))(t) = Lambda[v]-(beta[v]*i[h](t)+mu[v])*s[v](t), (D(i[v]))(t) = beta[v]*i[h](t)*s[v](t)-mu[v]*i[v](t);
##
awal := s[h](0) = 100, i[h](0) = 0, r[h](0) = 0, m[h](0) = 0, s[v](0) = 1000, i[v](0) = 0;
##
sol := dsolve({awal, sist}, numeric, maxfun = 10^6);
##
odeplot(sol, [t, s[h](t)], 0 .. 3*10^5);
## The next graph shows that i[h] is the zero solution:
odeplot(sol, [t, i[h](t)], 0 .. 100,thickness=3); #Dead zero
## The only other nonzero solution:
plots:-odeplot(sol, [t, s[v](t)], 0 .. 10^2,thickness=3);
## From the odes and the mostly homogeneous initial conditions it follows that
##  all but s[h] and s[v] must be identically zero:
##
for j from 1 to 6 do sist[j] end do;
awal;

Kitonum must be using the default setting of kernelopts(floatPi=true) in Maple 2017.
That means that Pi becomes its decimal approximation when in contact with a floating point number as in
6.557377048*Pi;
That kernelopts setting is not available in Maple 17 or 18.
###
It is available as an option in Maple 2015, but undocumented and the default is false.
In Maple 2016 and 2017 it is documented and the default is true.
I have set kernelopts(floatPi=false) in my maple.ini file. That is why I noticed: I got the same as you in my Maple 2017.
###
Another comment: When using inequalities in assume it is deduced that the quantity is assumed real.

In Maple 17 you can do:

restart;
assume(d>0):
assume(-0.01 < a, a < 0):
sys:={-800*Pi*a*cos(6.557377048*Pi*(3.470797713+d))/(a+1)^3 = -.9396060697, 800*Pi*a*sin(6.557377048*Pi*(3.470797713+d))/(a+1)^3 = -.3238482794};
## Extra step:
sys:=evalf(sys);
solve(sys, {a,d}, useassumptions=true, AllSolutions=true);

1. Put semicolon after end if. There are 3 of those, only for the first 2 of the 3 is the semicolon important though.
2. You need eval(..., list of equations). You have eval( ..., list of numbers).
    So use eval(XXX, [T,E,W,p]=~op_value) in two places.
3. Hessian should be VectorCalculus:-Hessian.
There is at least one more problem although the loops run for a while.
They end with the error: Array index out of range.

You should contact Customer Service directly:
http://www.maplesoft.com/support/

One simple advice I gave my students when separating variables in an ivp was to determine the arbitrary constant as early as possible. After observing that y(x)=0 is not interesting here, we would do:
dy/dx = x^(1/2)*y^(1/2)
y^(-1/2)*dy = x^(1/2)*dx
y^(1/2) = 1/3*x^(3/2) + C
Now determine C from y(0)=1 to be C = 1.
Squaring gives us the solution thereby avoiding the solve problem later.
It turns out that giving an optional argument usesolutions = "particular" or "particular via integrating factors"  singles out the correct solution:

sol := dsolve({ode, y(0) = 1},usesolutions="particular via integrating factors");
sol := dsolve({ode, y(0) = 1},usesolutions="particular");
## Trying all the possibilities mentioned in the help:
L:=["particular", "particular via integrating factors", "particular via symmetries", "particular via symmetries and general", "general", "general and particular", "general and particular via symmetries", "general and particular via integrating factors"];
#
for S in L do printf("Using %s\n",S); dsolve({ode, y(0) = 1},usesolutions=S) end do;

Notice that the first solution from

sol := dsolve({ode, y(0) = 1});
is in fact a solution of the ode, but not of the initial value problem.
It is a solution on the interval x = 3^(2/3) .. infinity:

simplify(odetest(sol[1],ode)) assuming x>3^(2/3);

I did your job in 1D input (aka Maple input):
 

restart;
varphi := x-> (m*omega/(Pi*`ℏ`))^(1/4)*exp(-(1/2)*m*omega*x^2/`ℏ`) ; 
solve(1 = abs(A)^2*(int(abs(varphi(x))^2, x = -infinity .. infinity)), A);
simplify([%]) assuming positive;

Answer: [1,-1]

Many errors can be caught by using a try statement as ThU is saying.
I give two examples below, where one is caught (the division by zero in the procedure p), but the other is not (the recursive assignment in q). My attention was called to the latter recursive example just recently, but my belief is that this try statement should work in the situation you mention.
 

restart;
p:=proc(x) 1/x end proc;
for i from -3 to 3 do
  try
    res:=p(i)
  catch:
    printf("Error at i = %d\n",i)
  end try;
  res
end do;
## Second example:
q:=proc(x) global E; E:='E'; if x=E then x:=x+1 else x end if end proc;
L:=[A,B,C,D,E,F,G]; #All unassigned!
res:='res':
for i in L do
  try
    res:=q(i)
  catch:
    printf("Error at i = %a\n",i)
  end try;
  res
end do;
## We see that the loop ends with D, so the uncatchable error occurred with E (not surprisingly considering the special status it has in q).
## We can, however, continue "manually" in the same session starting with F:
res:='res':
for i in L[-2..-1] do
  try
    res:=q(i)
  catch:
    printf("Error at i = %a\n",i)
  end try;
  res
end do;

If E is assigned a numeric value before the loop (e.g.  E:=5) then there won't be any errors, but I wanted to provoke a disastrous error.

Below I used the same initial condition as Kitonum.
It should be pointed out that method=classical by default is Euler's method (forward Euler).

restart;
res:=dsolve({diff(y(x),x) = -2*x*y(x) + 1, y(0)=2});
resN:=dsolve({diff(y(x),x) = -2*x*y(x) + 1, y(0)=2},numeric);#default method=rkf45
sol:=h->dsolve({diff(y(x),x) = -2*x*y(x) + 1, y(0)=2}, numeric, method=classical, stepsize=h);
plots:-odeplot(sol(.1),[x,y(x)],0..1); p1:=%:
plots:-odeplot(sol(.05),[x,y(x)],0..1,color=blue); p2:=%:
plot(rhs(res),x=0..1,color=green,thickness=3); p3:=%:
p4:=plots:-odeplot(resN,[x,y(x)],0..1,color=black); p4:=%:
plots:-display(p1,p2,p3,p4);

Notice that the exact solution (thick green) and the rkf45 numerical solution (black) are indistinguishable. That is my reason for having thickness=3 for the green one; otherwise that wouldn't be seen.

I had success when replacing bcs4 and bcs5 by
bcs4 := F1(inf) = 0, F3(inf) = 0, F5(inf) = 0;
bcs5 := G1(inf) = 0, G3(inf) = 0, G5(inf) = 0;
and with inf:=8.
I used an approximate solution as an extra argument to dsolve:
 

approxsoln = [F1(eta) = eta*(inf-eta), F3(eta) = -eta*(inf-eta), F5(eta) = -eta*(inf-eta), G1(eta) = exp(-eta), G3(eta) = eta*(inf-eta), G5(eta) = eta*(inf-eta), H1(eta) = -tanh(eta), H3(eta) = tanh(eta), H5(eta) = tanh(eta)]

Using inf:=10 I had success when I additionally used initmesh=128 as in

inf:=10;
R := dsolve({Eq1, Eq2, Eq3, Eq4, Eq5, Eq6, Eq7, Eq8, Eq9, bcs1, bcs2, bcs3, bcs4, bcs5}, numeric, output = listprocedure, approxsoln = [F1(eta) = eta*(inf-eta), F3(eta) = -eta*(inf-eta), F5(eta) = -eta*(inf-eta), G1(eta) = exp(-eta), G3(eta) = eta*(inf-eta), G5(eta) = eta*(inf-eta), H1(eta) = -tanh(eta), H3(eta) = tanh(eta), H5(eta) = tanh(eta)], initmesh = 128);

Supplying an approximate solution is a good idea in this case where the one found by dsolve itself is too simple. It uses
[F1(eta) = 0, F3(eta) = 0, F5(eta) = 0, G1(eta) = -(1/inf)*eta+1, G3(eta) = 0, G5(eta) = 0, H1(eta) = 0, H3(eta) = 0, H5(eta) = 0]