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]

With a and b as in your worksheet the simplest way seems to be:

Equate(a.b,<(0$4)>);
res:=solve(%);
eval(a.b,res); #Just a check

But you probably meant to require that the unknowns are real. Correct?

If so then you will have 8 linear equations with 7 unknowns. That system is likely not to have any solution. Indeed it appears so:

eqs:=Equate(a.b,<(0$4)>);
eqs1:=evalc(Re~(eqs));
eqs2:=evalc(Im~(eqs));
solve({op(eqs1),op(eqs2)});

Since your numbers are decimal numbers (floats) the existence of a solution to your overdetermined system will be sensitive to roundoff in input etc.

Thus you may want to settle for the best substitute for a solution in the least squares sense:

nms:=convert(indets(b),list);
A,B:=LinearAlgebra:-GenerateMatrix([op(eqs1),op(eqs2)],nms);
X:=LinearAlgebra:-LeastSquares(A,B);
LinearAlgebra:-Norm( A . X - B ,2); # good if small
eval(a.b,nms=~convert(X,list)); #Another check

You can solve numerically.
The heat equation immediately smoothes the initially given function and your boundary conditions bc:=f(t,0)=0,f(t,1)=1 together with the initial condition ic:=f(0,x)=piecewise(x=0,1,0) disagree at (t,x) = (0,0) and (t,x) = (0,1). Thus you can expect to see problems initially in the numerical solution.
I wonder, however, if you didn't intend ic to be ic:=f(0,x)=piecewise(x=1,1,0) instead? That wouldn't help for an exact solution, but you won't see any numerical problems:
 

sol:=pdsolve(pde,{bc,ic},numeric,spacestep=0.05,timestep=0.001);
sol:-animate(t=0.4,frames=100);

Use numerical solution:
 

restart;
ode:=diff(y(x),x,x)=5*exp(-10/diff(y(x),x)); 
bcs:=y(0)=0,y(15)=2;
res:=dsolve({ode,bcs},numeric);
plots:-odeplot(res,[x,y(x)]); # Looks like a straight line y=2*x/15
plots:-odeplot(res,[x,y(x)-2*x/15]); # basically 0
odetest(y(x)=2*x/15,ode); # Almost zero, but not quite

I tried the model suggested by Carl Love and modified by Kitonum in NonlinearFit in the Statistics package.
 

restart;
M:=ImportMatrix("C:/Users/Bruger/Downloads/Book1.xlsx",datatype=float[8]);
M[1..10,..];#Having a look
mx:=max(M[..,1]); #maximal x-value
F:=Statistics:-NonlinearFit(A*x^n*sin(B*x),M,x,initialvalues=[A=0.015,B=280,n=1.5],output=solutionmodule);
F:-Results();
res:=F:-Results("leastsquaresfunction");
sq1:=F:-Results("residualsumofsquares"); #For comparison with Kitonum's function
f:=x->0.015*x^1.5*sin(280*x); #Kitonum
Y:=f~(M[..,1]);
sq2:=add((Y[i]-M[i,2])^2,i=1..numelems(Y)); #"residualsumofsquares" for f.
p1:=plot(M,style=point):
p2:=plot(res,x=0..mx,color=blue,numpoints=30000):
plots:-display(p1,p2,size=[1800,default]);
p3:=plot(f(x), x=0..mx,color=blue,numpoints=30000):
plots:-display(p1,p3,size=[1800,default]);
plots:-display(p1,p3,size=[1800,default],view=[mx-.1..mx,default]);
plots:-display(p1,p2,size=[1800,default],view=[mx-.1..mx,default]);
################# Introducing 2 extra parameters:
F2:=Statistics:-NonlinearFit(A*x^n*sin(B*x^m+phi),M,x,initialvalues=[A=0.015,B=280,n=1.5,m=1,phi=0],output=solutionmodule);
res2:=F2:-Results("leastsquaresfunction");
F2:-Results("residualsumofsquares"); #0.0008
p22:=plot(res,x=0..mx,color=blue,numpoints=30000):
plots:-display(p1,p22,size=[1800,default]);

The result 'res' was

Kitonum's function plot in some way looks better by hitting the tops, but the residual sum of squares is 10 times larger than that from NonlinearFit. Also the closeups at the end:

plots:-display(p1,p2,size=[800,default],view=[mx-.1..mx,default],caption="NonlinearFit");
plots:-display(p1,p3,size=[800,default],view=[mx-.1..mx,default],caption="Kitonum's function");

are revealing: