The short answer is no! The signs of b(x) and a(x) are the ones that are important, not x.

If you look at the help page for arctan you will find this statement:
"For real arguments x, y, the two-argument function arctan(y, x), computes the principal value of the argument of the complex number  x + I y, so  -Pi < arctan(y, x) and arctan(y, x) <= Pi."
 You could try this:
 

arctan(y,x) assuming x>0, y::real;
arctan(y,x) assuming x<0, y>=0;  
arctan(y,x) assuming x<0, y<0;
arctan(y,0) assuming y>0;
arctan(y,0) assuming y<0;
arctan(0,0);

With the default Digits=10 I actually get an error:

Error, (in Optimization:-LPSolve) invalid input parameter to external routine

Raise Digits to 25 and you will get the result you expect. So do
Digits:=25;
before using Maximize.

In doing
evalf( exp(1) );
you are getting a decimal number approximation to exp(1). With Digits at 10 ( the default) that number is 2.718281828 .
You may try
evalf[50](exp(1));
and you get 2.7182818284590452353602874713526624977572470937000 .
Try replacing 50 by e.g. 500.
Maple will not give you an approximation of a number like exp(1) or Pi unless you ask for it either directly as above or indirectly as in
exp(1.);
exp(1.0); #same thing
where I have replaced the integer 1 with the decimal number 1.0 (or just 1. ).
I^2 evaluates to -1 exactly.
But try
I^2.0;
sqrt(2)^2 evaluates to 2 exactly.
But try
sqrt(2.0)^2;
Don't use decimal numbers if you are not satisfied with an approximation.

If you don't need the various statistical information that can be supplied by LinearFit you can use LSSolve:
 

restart;
pts:=[[0,0],[1,13],[2,21],[6,45],[12,54],[15,77]];
X:=evalindets[2](pts,list(realcons),1,op);
Y:=evalindets[2](pts,list(realcons),2,op);
f:=x->a*(x-2)+21; # y-21 = a*(x-2)
d:=f~(X)-Y;
Optimization:-LSSolve(d);

The successful methods are ftoc and ftocms (ftoc stands for The Fundamental Theorem of Calculus).
If you find an antiderivative:

F:=int(sin(K*x)*sin(L*x)*cos(M*x), x);

you will see that
F := (1/4)*sin((K-L-M)*x)/(K-L-M)+(1/4)*sin((K-L+M)*x)/(K-L+M)-(1/4)*sin((K+L-M)*x)/(K+L-M)-(1/4)*sin((K+L+M)*x)/(K+L+M)

thus this is only valid if K-L-M, K-L+M, K+L-M, and K+L+M are all different from zero.
In your case K-L-M = 14-2-12 = 0, so the generic F cannot be used.
### WORKAROUND:
You can use that if a parameter is present in the limits then the option AllSolutions will give all cases.
 

restart;
A:=Int(sin(K*x)*sin(L*x)*cos(M*x), x = 0 .. Pi);
## We may assume that K<>0 since otherwise A = 0 no matter the values of L and M.
B:=IntegrationTools:-Change(A,x=t/K,t); # Assuming that K<>0
f,r:=op(op(2,B));
int(f,r,AllSolutions)/K assuming integer;
C:=simplify(%);
eval(C,{K=14,L=2,M=12});

Your system may have no solution or at least such a solution is very difficult to find for the following reason.
Solving for the highest derivatives gives you a system whose right hand sides has a denominator that will take the value 0 in the open interval (0, L). That just follows from the boundary conditions as you will see here:

restart;
Eq1 := diff(F(eta), eta, eta, eta)-phi1*((diff(F(eta), eta))^2-(diff(F(eta), eta)+diff(G(eta), eta))*(diff(F(eta), eta, eta)))+alpha*(1-phi)^2.5*((diff(F(eta), eta, eta))^2+2*(diff(F(eta), eta))*(diff(F(eta), eta, eta))-(diff(F(eta), eta)+diff(G(eta), eta))*(diff(F(eta), eta, eta, eta, eta))+A*(2*(diff(F(eta), eta, eta, eta))+(1/2)*eta*(diff(F(eta), eta, eta, eta, eta))))-A*phi1*(diff(F(eta), eta)+(1/2)*eta*(diff(F(eta), eta, eta))) = 0; Eq2 := diff(G(eta), eta, eta, eta)-phi1*((diff(G(eta), eta))^2-(diff(F(eta), eta)+diff(G(eta), eta))*(diff(G(eta), eta, eta)))+alpha*(1-phi)^2.5*((diff(G(eta), eta, eta))^2+2*(diff(G(eta), eta))*(diff(G(eta), eta, eta))-(diff(F(eta), eta)+diff(G(eta), eta))*(diff(G(eta), eta, eta, eta, eta))+A*(2*(diff(G(eta), eta, eta, eta))+(1/2)*eta*(diff(G(eta), eta, eta, eta, eta))))-A*phi1*(diff(F(eta), eta)+(1/2)*eta*(diff(F(eta), eta, eta))) = 0;
###
bcs:=F(0)=0,D(F)(0)=1,G(0)=0,D(G)(0)=p, D(F)(L)=0, D(G)(L)=0, F(L)=0, G(L)=0;
###
phi1 := 1.2; alpha := 2; A := 1.5; phi:=0.1; p:=0.1; L:=8;
sys:=solve({Eq1,Eq2},diff~({F,G}(eta),eta$4));
(denom@rhs)~(sys); # The same denominator in both equations
DN:=op(%);
eval[recurse](convert(DN,D),{eta=0,bcs}); # Value 1.690553637*10^9 > 0
eval[recurse](convert(DN,D),{eta=L,bcs}); # Value -9.221201656*10^9 < 0

Even if you change the late coming two boundary conditions won't alter the result of the last two lines: Try using

bcs6:=F(0)=0,D(F)(0)=1,G(0)=0,D(G)(0)=p, D(F)(L)=0,D(G)(L)=0;

instead of bcs in the last two lines above.

Try this:
 

restart;
f:=x->sqrt(2.0)*x;
y := sqrt(3.0); #Digits = 10
for n from 10 to 3 by -1 do
  [evalf[n](evalf[n](sqrt(2.0))*evalf[n](y)),evalf[n](f(y))]
end do;

The two numbers for all n = 3..10 are exactly the same.
And by the way, so are these:

for n from 10 to 3 by -1 do
  [evalf[n](sqrt(2.0)*y),evalf[n](f(y))]
end do;

Go to Tools/Options/Display and uncheck Show equation labels.
Apply Globally or Apply to Session.

restart;
f:=x->x^2-3;
x:=1;
for i from 2 to 5 by 1 do
    x:=evalf(x-f(x)/D(f)(x));
    print(x);
end do;
evalf(sqrt(3));

Besides what Kitonum has already pointed out, you must leave 'continuation' out. There is no continuation parameter in your system.

Furthermore, it appears that the space between Array and the list of numbers is being interpreted as multiplication (as it will be in 2D math input):
It should be:
 

output = Array([0., 0.5e-1, .10, .15, .20, .25, .30, .35, .40, .45, .50, .55, .60, .65, .70, .75, .80, .85, .90, .95, 1.00])

Finally, for plotting you only need to do

with(plots); 
odeplot(r1, [x, f(x)]);

No need for display (you missed an assignment there too, incidentally).
 

The problem is that when the procedure is called as in JD(2018, 1, 11) , then the first that happens is that all occurrences of YY, MM, and DD in the procedure body are replaced with the numbers 2018, 1, and 11 respectively.
The assignments YY := YY-1; MM := MM+12 thereby becomes
2018:=2018-1; 1:=1+12;
But needless to say, you cannot assign to numbers like 2018 or 1.
So change to
 

JD := proc (Y, M, DD)
  local A, B, YY, MM;
  YY:=Y; 
  MM:=M;
  if MM = 1 or MM = 2 then
    YY := YY-1; 
    MM := MM+12
  end if;
  A := floor((1/100)*YY);
  B := 2-A+floor((1/4)*A);
  evalf(floor(365.25*YY)+floor(30.6001*(MM+1))+DD+1720994.5) 
end proc;

Whether your procedure works as intended is another matter.

Here is an example:
 

restart;
ode:=diff(y(x),x)=y(x);
h:=0.1: #Stepsize
resRKF45:=dsolve({ode,y(0)=1},numeric,output=Array([seq(h*i,i=0..20)]));
resRK4:=dsolve({ode,y(0)=1},numeric,method=classical[rk4],stepsize=h, output=Array([seq(h*i,i=0..20)]));
A45:=resRKF45[2,1];
A4:=resRK4[2,1];
X:=A4[..,1];
Yd:=A4[..,2]-A45[..,2];
plot(X,Yd,caption=typeset("The difference between RK4 and RKF45\n with RK4 stepsize = ",h));

In this case you can also compare the two results with the exact solution y(x) = exp(x):
 

plots:-odeplot(resRKF45,[x,y(x)-exp(x)],0..2);
plots:-odeplot(resRK4,[x,y(x)-exp(x)],0..2);

@mrbayat There doesn't seem to be a difference between the mw and the mws file. Both have 2D-input, thus the mws-file isn't an mws-file.

But apart from that:
I executed your whole worksheet by using the !!! icon. Then I did:
 

indets(dsys,function);

What I saw includes weird constructions like X2(X2), diff(X2(X2), X2).

This must be cleared up.
The best is to remove X2 in N(X2) from the first line after restart so that you get:
 

Ws := (1/2)*N*k*T*(I1-3-2*log(J));

It appears as if a multiplication sign is missing in the line defining Wm. In any case it should be:

Wm := k*T*((J-1)*log(1-1/J)+chi-chi/J)/nu;

Finally, at the end you now get an error message saying that 7 boundary conditions are expected, but only two were provided. That is due to the fact that there are several (5) undefined constants in your system. Try

indets(dsys,name);

You get {L, T, X2, chi, k, nu}. So L, T, chi, k, nu needs to be given numerical values.

Have a look at PDEtools:-dchange.
For the present problem you can do as follows, where I have included initial conditions at t = t0.
First I'm doing ode and ics separately (more straightforward):
 

restart;
ode:=m*diff(x(t),t,t)+k*x(t)=0;
ics:=x(t0)=x0, D(x)(t0)=x1;
## Doing ode and ics separately:
PDEtools:-dchange({t=T*s, x(t)=X(s)*L},ode,[s,X]);
solve(%,{diff(X(s),s,s)});
subs(T^2=m/k,%);
ICS:=PDEtools:-dchange({t=T*s, x(t)=X(s)*L},{ics},[s,X]);
indets(ICS,function);
icsX:=subs(T=sqrt(m/k),solve(ICS,%));
##############
## Doing ode and ics at the same time:
res:=PDEtools:-dchange({t=T*s, x(t)=X(s)*L},{ode,ics},[s,X]);
## Separating the DE from the initial conditions
res1,res2:=selectremove(has,res,diff);
solve(res1,{diff(X(s),s,s)});
odeX:=op(subs(T^2=m/k,%));
indets(res2,function);
icsX:=subs(T=sqrt(m/k),solve(res2,%));

You may want to introduce s0 = t0*sqrt(k/m), X0 = x0/L, X1 = (x1/L)*sqrt(m/k).
While solving the new DE for the derivative(s) (here diff(X(s),s,s) ) is not necessary, it makes it easier to see (in general) which nondimensionalization to choose.
## You have introduced the constant L, which you can still choose to fit the situation. If x0 is never zero then L could be chosen as x0 so that X(s0) = X0 = 1.

Your second pade expression is when copied by me this:
 

pade*(theta, eta, [3, 3])

Obviously it should be:
 

pade(theta, eta, [3, 3])

With that change fsolve returns {b = -.26165340, d = .14500436}.