The exponential function exp, should be used like this exp(x), NOT like this e^(x).
I corrected that below.
You assign pi to 3.143. But you also use Pi. Below I have assumed that pi was intended as a (bad) approximation to Pi, so I have used Pi in all places.
So we have:
 

restart;
Digits:=50:
G := 6.6743*10^(-8);
a := 1.9501*10^24;
b := .3306;
c := 2.99792458*10^10;
d := 2.035;
#pi := 3.143; #Pi also appears below. I changed pi to Pi.
eps := 3.8220*10^35;
## g(r) and j(r) are not used in the following:
#g(r) = 1-s(r)/0.06123;
#j(r) = exp(-(1/2)*w(r))*(1-2*G*v(r)/(r*c^2))^(1/2); #Must use exp( ) not e^( )
sys1 := diff(v(r), r) = 4*Pi *r^2*eps/c^2, v(0)=0;
res1:=dsolve({sys1});
sys2:=diff(u(r), r) = -G*(eps+u(r))*(v(r)+4*Pi*r^3*u(r)/c^2)/(c^2*(r^2-2*G*r*v(r)/c^2)),u(0)=1.3668*10^34;
sys3:=diff(w(r), r) = 1.485232054*10^(-28)*(v(r)+4.450600224*10^(-21)*Pi *r^3*u(r))/(r^2-2*G*r*v(r)/c^2),w(0)=0;
## Now solve numerically
res2:=dsolve(eval({sys2,sys3},res1),numeric);
plots:-odeplot(res2,[[r,u(r)],[r,10^35*w(r)]],0..700000);
res2(688240);

The result of the last command is that w(r) = 0.0704472 ... and not as you "want"  w(688240)=-2.05684.
You cannot impose both conditions w(0)=0 and w(688240)=-2.05684.
 

You have a variety of choices. Here I give 3 methods:

restart; 
with(LinearAlgebra): 
A := Matrix(2,2,[-4,-2,3,1]); 
X:=<x1(t),x2(t)>;
ode:=diff~(X,t)=A.X; #The ode 
##Finding eigenvalues and corresponding eigenvectors (the columns in P): 
Lambda,P:=Eigenvectors(A); ## The general solution is
X = add(c[i]*exp(Lambda[i]*t)*Column(P,i),i=1..2); #FIRST
## That (in another form) could have been found by dsolve: 
dsolve(ode); 
X=subs(%,X); #SECOND
## You also have available MatrixExponential:
E:=MatrixExponential(A,t);
## The general solution in a third form: 
X = E.<c1,c2>; #THIRD

If you are not using a recent version of Maple, then dsolve won't take vector input.
In that case, replace ode by

ode1:=Equate(op(ode));
## and then as above:
dsolve(ode1);
X=subs(%,X);

Since it appears that you want to compare the numerical solution with the exact one you can try this:

restart;
Digits:=15;
ode := diff(w(x), x$6)+diff(w(x), x$4)+diff(w(x), x, x)-w(x) = -90;
bcs:=w(0) = 0, w(1) = 0, (D@@2)(w)(0) = 0, (D@@2)(w)(1) = 0, D(w)(0) = 0, D(w)(1) = 0; 
res:= dsolve({ode,bcs}, numeric);
plots:-odeplot(res,[x,w(x)]);
pol:=DEtools[de2diffop](lhs(ode),w(x),[R,x]); #The characteristic polynomial
rts:=solve(pol=0,R); #Roots not "nice"
evalf(rts);
## An impatient attempt at finding the exact solution:
sol:=timelimit(30,dsolve({ode,bcs})); #30 seconds expires
## I gave up on that also because the solution would take up a lot of space.
####
solG:=dsolve(ode): #The general solution
## Now imposing the boundary conditions:
eq1:=eval(rhs(solG),x=1):
eq2:=eval(diff(rhs(solG),x),x=1):
eq3:=eval(diff(rhs(solG),x,x),x=1):
eq4:=eval(rhs(solG),x=0):
eq5:=eval(diff(rhs(solG),x),x=0):
eq6:=eval(diff(rhs(solG),x,x),x=0):
## Use fsolve; it is fast:
S:=fsolve({eq1,eq2,eq3,eq4,eq5,eq6});
S:=simplify(fnormal(S,10));
sol:=simplify(evalc(evalf(eval(solG,S))));
plot(rhs(sol),x=0..1);
plots:-odeplot(res,[x,w(x)-rhs(sol)]); #Comparison

The reason is that there is no solution. You may also try some linear algebra:

M,b:=LinearAlgebra:-GenerateMatrix([A1, A2, A3, A4, A5, A6],[C1, C2, C3, C4, C5, C6]);
LinearAlgebra:-LinearSolve(M,b); #Inconsistent system, so no solution
## Details:
LinearAlgebra:-GaussianElimination(<M|b>);
### QUESTION:
What does the second part have to do with the first?
Consider
ode:=diff(w(x), x, x, x, x, x, x)+diff(w(x), x, x, x, x)+diff(w(x), x, x)+(1-2)*w(x) = -90;
dsolve(ode);
evalf(%);
## So why are you looking at
w(x)= C1*sinh(x)+C2*cosh(x)+C3*sin(x)+C4*cos(x)+C5*sin(x)+C6*cos(x)+90;
?
Notice that regardless of that question, the appearance of sin(x) and cos(x) twice each is strange!

Begin by deleting maxfun as the error message says.
Then you could try using approxsoln=[psi(u)=...] as an option.

My own experiments indicate that a solution looks like this:

Use the view option in Maple 2015 (which you are using, it is not needed in Maple 2016):

 plot3d(2*x/(x^2+y^2), x = -10 .. 10, y = -10 .. 10,view=-5..5);

For commenting out several lines use   (*  ....   *)  as in

(*

line 1
line 2
line 3
*)
 

You have a system of two odes that is third order in f and first order in theta.

Since 3+1 = 4, you need 4 boundary conditions.

In your other version, which I deleted, the problem is the same.

In the help page in Maple 2016, dsolve/details, you find under Optional Arguments:

'parametric'
To request the use of only the parametric solving scheme when solving a single first order ODE, possibly of high degree in dy/dx. Note however that by default dsolve tries to remove the parameter used during the solving process. To keep the parameter, specify the optional argument 'implicit' together with 'parametric'.

In the Examples section on that page you find an example (look for ode[4]);  you could also take a look at the link given there in the text.

Changing your input from vector to list makes your corrected code work:

restart;
z:=A*[x,y];
A:=1.0;
plots:-fieldplot(z,x=0..1,y=0..1);

Notice that all examples given in the help page for fieldplot has a list as the first argument. Notice also that the help page states:

A typical call to the fieldplot function is fieldplot(f, x=a..b, y=c..d), where f is a list of two expressions in x and y. a..b and c..d specifies the horizontal and vertical ranges of the field. The result is a two-dimensional vector field with the vector evaluated at (x, y) located at this point.
(emphasis added).

It is also mentioned that a VectorCalculus:-Vectorfield can be used, and indeed it works too:

restart;
z := VectorCalculus:-VectorField( A*<x, y>, 'cartesian'[x,y] );
A:=1.0;
plots:-fieldplot(z,x=0..1,y=0..1);
## Also this version works fine:

restart;
VectorCalculus:-BasisFormat(false);
z := VectorCalculus:-VectorField( A*<x, y>, 'cartesian'[x,y] );
A:=1.0;
plots:-fieldplot(z,x=0..1,y=0..1);

Try this:
First case:
 

restart;
ode:=diff(y(x),x,x)+diff(y(x),x)+y(x)=F;
bcs:=y(0)=1,y(2)=1;
F  := piecewise(x<=1, -x, x>1, 2*x+(x-1)^2);
sol:=dsolve({ode,bcs}); #That is not your solution
plot(rhs(sol),x=0..2); pE:=%:
plot(diff(rhs(sol),x),x=0..2);
plot(diff(rhs(sol),x,x),x=0..2,discont);
As:=dsolve({ode, bcs} ,numeric, maxmesh=500,abserr=1e-3); 
plots:-odeplot(As,[x,y(x)-rhs(sol)]); #Difference between numerical and exact solution: OK considering abserr

Now continue with the second case:

####
F := piecewise(x<=1, x^2+x+2, x>1, -x^2+x);
ode;
sol:=dsolve({ode,bcs});
plot(rhs(sol),x=0..2); 
plot(diff(rhs(sol),x),x=0..2);
plot(diff(rhs(sol),x,x),x=0..2,discont);
As:=dsolve({ode, bcs} ,numeric, maxmesh=4000,abserr=1e-3); 
plots:-odeplot(As,[x,y(x)-rhs(sol)]);

In the first case your exact "solution" on the interval 0..2 is less smooth than the solution provided by Maple's dsolve. Notice that if you rewrite the 2nd order equation as a first order system, then in your case one of the dependent variables (the one that was y'(x)) is not continuous. It should definitely be a requirement for a solution to a first order system on an interval a..b that the functions are (at least)  continuous on the whole interval a..b.
### To show the reason for this natural requirement, I shall show that if you give up continuity of the first derivative in your first case, then you have infinitely many solutions of your bvp problem:
 

F:= piecewise(x<=1, -x, x>1, 2*x+(x-1)^2);
sol1:=dsolve({ode,y(0)=1}) assuming x<1;
sol1:=subs(_C2=C1,sol1);
sol2:=dsolve({ode,y(2)=1}) assuming x>1;
sol2:=subs(_C2=C2,sol2);
eval(rhs(sol1),x=1)=eval(rhs(sol2),x=1); #One equation, two unknowns!
##So you can pick C2 to be anything, then C1 is determined.
##Solving for C1 in terms of C2:
solve(%,{C1});
sol11:=eval(sol1,%);
solT:=y(x)=piecewise(x<=1,rhs(sol11),rhs(sol2)); # "solution" on 0..2
plots:-animate(plot,[rhs(solT),x=0..2],C2=-1..1);
plots:-animate(plot,[diff(rhs(solT),x),x=0..2,discont],C2=-1..1);

The animation of y(x) above:

I shall assume that all variables are positive. You can do with less, but try:

Ex3 := Pi*sqrt((4/mu+4/r)*exp(-r/mu));
simplify(Ex3) assuming positive;
combine(%) assuming positive;

You have eq3= ... and eq4= ... .
You need assignment, so use eq3:= ... and eq4:= ....

I tried dsolve/numeric:

res:=dsolve(ODE union BC1 union {f(0) = 0, fp(0) = .2*fpp(0),h(0) = 0,i(0) = 0,g(0) = 1+.2*gp(0)},numeric);
Why do you want to use a shooting method?
Anyway, if you do, try this:

1. replace ODE by ODE1 given by
ODE1:=solve(ODE,indets(ODE,specfunc(diff)));

2. Replace IC by IC2 given by the explicit initial conditions:
IC2 := {f(0) = 0, fp(0) = gamma1*delta, fpp(0)=delta, g(0) = 1+gamma1*beta, gp(0) = beta, h(0) = 0, hp(0) = beta1, i(0) = 0, ip(0) = beta2, fppp(0) = alpha};

3. Then do:

S := shoot(ODE1, IC2, BC1, FNS, [alpha = .1, beta = .2, beta1 = .3, beta2 = .4,delta=0]);

I tried to insert a print statement in both obj1 and obj2, here only shown in obj1:
 

obj1:=proc(u1,u2)
local z1,s1,s2,z2,res;
global sol1;
sol1('parameters'=[1,0,u1]);
z1 := sol1(.5);
s1:=subs(z1,y1(t));
s2:=subs(z1,y2(t));
sol1('parameters'=[s1,s2,u2]):
z2:=sol1(0.5);
res:=-subs(z2,y2(t));
print(res);
res
end proc;

Then I tried your worksheet. That showed that fdiff in computing Hess1 uses only 18 digits (i.e. the global setting of Digits + 3), whereas in Hess2 twice as many digits are used.Then I tried a third version obj3, where sol1 is made local. That mostly defeats the purpose of using the parameters option in dsolve, I suppose, but here it is:

obj3:=proc(u1,u2)
local z1,s1,s2,z2,res,sol1;
sol1 := dsolve({sys, y1(0) = alpha, y2(0) = beta}, type = numeric, 
'parameters' = [alpha, beta, u],maxfun=0,range=0..0.5):
sol1('parameters'=[1,0,u1]);
z1 := sol1(.5);
s1:=subs(z1,y1(t));
s2:=subs(z1,y2(t));
sol1('parameters'=[s1,s2,u2]):
z2:=sol1(0.5);
res:=-subs(z2,y2(t));
print(res);
res
end proc;

I then got (almost) exactly Hess2, when doing:

Hess3:=Matrix(2,2):
for i from 1 to 2 do 
  for j from 1 to 2 do 
    Hess3[i,j]:=fdiff(obj3,[i,j],u0,workprec=1.0)
  od
od:
Hess3;

## I tried Digits:=10, then I got Hess1 roughly equal to Hess2.