#Answer edited to raise Digits initially inspired by comparison with the numerical solution.
Maple can solve the system without the boundary conditions as you also observed. Then you can procede with determining the 6 arbitrary constants yourself using solve (or fsolve).
## Start the session with Digits:=25 (something much higher than 10).
restart;
Digits:=25:
..........
..........
..........
res:=dsolve({sys1}):
indets(res,name);
indets(res,specfunc(anything,Int)):
nops(%); #2 integrals Int
res2:=value(res);
indets(res2,specfunc(anything,{Int,int})); #Turning Int into int evaluates the integrals successfully
bcs2 := X11(0) = 0, X11(L) = 0, X22(0) = 0, X22(L) = 0, X33(0) = 0, X33(L) = 0;
subs(x=0,bcs2,res2): eq1:=%;
subs(x=10,bcs2,res2): eq2:=%;
evalf(eq1 union eq2);
solve(%,{_C1,_C2,_C3,_C4,_C5,_C6});
sol:=eval(res2,%);
plot(subs(sol,[10000*X11(x),X22(x),X33(x)]),x=0..10);
#Notice that I have multiplied X11 by 10000

To see the definition of Ci and Si just do ?Ci in Maple

A float is a decimal number like 1.23456.

But your problem is simply that the integral v2 is divergent.

Try this one, which has the same problem at 0:

int((1-sin(y*u))/u^3, u = 0 .. 1);
The output is infinity.

In your worksheet the output at the bottom (where y=1) is


showing the same thing: the integral having u^3 in the denominator is divergent.

If as you say the other people's software gives an answer that is not infinity then they are either considering another problem or they should buy Maple.

#################
I am beginning to think that you probably are missing a parenthesis, so that you actually wanted

f:=((1-exp((-t*u^2)/(1+alpha*u^2)))*sin(y*u))/(u^3);

and not the one you got. If so there is no divergence and the other people probably have the correct f.

You have theta(3)=0 as boundary condition and h=0.3. Number of steps = 10.
Thus you need theta[10]:=0 NOT theta[3]:=0.

Use this for d6 instead of the one you got:
d6 := (theta[i+1]-2*theta[i]+theta[i-1])/h^2

In the loop assigning to the theta equations the loop parameter has to be i since that index was used in ode2. So change the loop to:
for i to 9 do eq[i] := ode2 end do;

Inspired by the results of dsolve, which shows that theta decreases from 1 to 0, better guess values would be as in:

soln2 := fsolve({eq[1], eq[2], eq[3], eq[4], eq[5], eq[6], eq[7], eq[8], eq[9]}, {theta[1] = 1, theta[2] = 1, theta[3] = 0, theta[4] = 0, theta[5] = 0, theta[6] = 0, theta[7] = 0, theta[8] = 0, theta[9] = 0});

You could also try without guess values.

You are presenting a strange problem. Initial and boundary conditions should be given along curves not just at a single point (x0,y0). So I tried looking at the image you linked to and did the following in order to recover the equations in the image.
There is a difference:
Some squares appearing in your equations don't appear in the image. Is that due to a bad image?

restart;
eq1:=diff(U2(x, y), y)*((L-U1(x, y))^2+(D2+tan(alpha)*(L+U1(x, y)))^2)+.5*U2(x, y)*(2*tan(alpha)*U1(x, y)-D2) = 0;
EQ1:=subs(U2(x,y)=diff(U1(x, y), y),U1=h,eq1);
#Except for some missing squares this is equivalent to one of the image equations
eq2:=diff(U3(x, y), x)*((L-U1(x, y))^2+D1^2*cos(alpha+2*U2(x, y))^2)+.5*U3(x, y)*D1*cos(alpha+2*U2(x, y))+.5*D1*sin(alpha+2*U2(x, y))*(L-U1(x, y)) = 0;
EQ2:=subs(U3(x,y)=diff(U1(x, y), x),U2(x,y)=diff(U1(x, y), y),U1=h,eq2);
#Except for some missing squares this is equivalent to the other image equation
conds:={h(x0,y0)=K1,D[1](h)(x0,y0)=K2,D[2](h)(x0,y0)=K3};
#These conditions are strange.

##The following has been edited.
#We can consider EQ1 an ordinary differential equation in the variable y for each fixed x. If we only have the conditions conds EQ1 will determine h(x0,y) as a function of y. We cannot know anything about h(x,y) for other values of x.
The function h(x0,y) can be found like this, where we take x0=0:
We need to solve numerically so we must pick numeric constants.
If we pick both K1=0 and K3=0 we get the rather boring zero solution, so I picked K3=1.
EQ1d:=subs(h(x,y)=h(y),EQ1);
res:=dsolve({eval(EQ1d,{L=1,D2=1,alpha=1}),h(0)=0,D(h)(0)=1},numeric,output=listprocedure);
H,H1:=op(subs(res,[h(y),diff(h(y),y)]));
plots:-odeplot(res,[y,h(y)],y=0..500);
plots:-odeplot(res,[y,diff(h(y),y)],y=0..500);

You could determine h(x,y) from EQ1 if you had information like h(x,0)=f(x),D[2](h)(x,0)=g(x) where f and g are known functions. You just do as above for each x in a range. Using the parameters option in dsolve with parameters=[x] would do it easily.

But then what about EQ2?

I downloaded the file. Opened it in Maple 18.
Up came a statement.
march('open',"C:\\Users\\Preben\\AppData\\Local\\Temp\\CHEMISTRY.mla");
After accepting execution I wrote:
with(CHEMISTRY);
Then I could look at the exported procedures like in this example:

showstat(DeduceExperiments);

If you want to examine the locals you could do:
kernelopts(opaquemodules=false):

eval(CHEMISTRY:-GKBName);

showstat(CHEMISTRY:-FoundGoals);

The sequence of locals:

L:=op(3,eval(CHEMISTRY));

Checking which are procedures:
type~([L],procedure);
#or better
select(type,[L],procedure);
#Similarly
select(type,[L],table);
eval(L[3]);

Not that it really affects the comparison, but in the finite difference scheme you actually are considering the interval 0..5, not 0..10. Thus the code in P2.mw should have the interval 0..5, i.e the 10 in bcs should be 5.

When plotting sols in P1.mw:
plot(sols);

the curve doesn't reflect the fact that the first derivative at eta=5 should be zero.

The problem may be eq[11] and eq[12], see below.

I suppose that the left hand side of eq[12] is an adequate discrete approximation of the second derivative at the end eta=5 (or maybe in fact eta=4.5)? See revised worksheet at the bottom.

Notice that eta[10]=5. I suppose that f[10] corresponds to that eta.
Thus I tried changing eq[11] and eq[12]:

eq[11] := (f[10]-f[9])/h = 0;
eq[12] := (f[10]-4*f[9]+6*f[8]-4*f[7]+f[6])/h^2 = 0;

Now at least plot(sols) does reflect that f has derivative zero at eta=5.

There may be other problems, but the nagging question is: Is there actually more than one solution? And you just found another than dsolve?

####A revised worksheet P1.mw :

restart;
h := .5; k_1 := 0.9e-1;
ode := (diff(f(eta), eta))^2-f(eta)*(diff(f(eta), eta, eta)) = diff(f(eta), eta, eta, eta)+k_1*(2*(diff(f(eta), eta))*(diff(f(eta), eta, eta, eta))+(diff(f(eta), eta, eta))^2-f(eta)*(diff(f(eta), eta, eta, eta, eta)));
bcs := f(0) = 0, (D(f))(0) = 1, (D(f))(10) = 0, ((D@@2)(f))(10) = 0;
d1 := (f[i+1]-f[i-1])/(2*h);
d2 := (f[i+1]-2*f[i]+f[i-1])/h^2;
d3 := (f[i+2]-2*f[i+1]+2*f[i-1]-f[i-2])/(2*h^3);
d4 := (f[i+2]-4*f[i+1]+6*f[i]-4*f[i-1]+f[i-2])/h^4;
## You don't use these in fact!
#fd := (f[i+1]-f[i])/h;
#bd := (f[i]-4*f[i-1]+6*f[i-2]-4*f[i-3]+f[i-4])/h^4;
ode1 := d1^2-d2*f[i] = d3+k_1*(2*d1*d3+d2^2-d4*f[i]);
eta[0] := 0;
for j from 0 to 9 do eta[j+1] := eta[j]+h end do;
#Using boundary conditions:
f[0] := 0;
eq[10] := (f[1]-f[0])/h = 1;
eq[11] := (f[10]-f[9])/h = 0;
eq[12] := (f[11]-2*f[10]+f[9])/h^2 = 0; ##Notice change
#I changed from 0..9 to 1..9:
for i to 9 do eq[i] := ode1 end do:
indets({eq[1], eq[2], eq[3], eq[4], eq[5], eq[6], eq[7], eq[8], eq[9], eq[10], eq[11], eq[12]});
#With guessing values taken from dsolve in the other worksheet:
soln:=fsolve({eq[1], eq[2], eq[3], eq[4], eq[5], eq[6], eq[7], eq[8], eq[9], eq[10], eq[11], eq[12]}, {f[-1]=0, f[1]=0.4, f[2]=0.65, f[3]=0.8, f[4]=0.9, f[5]=0.96, f[6]=0.99, f[7]=1.01, f[8]=1.02, f[9]=1.028, f[10]=1.029, f[11]=1});
assign(soln);
sols := [seq([eta[i], f[i]], i = 0 .. 10)];
plot(sols);

You have three syntax problems:

1. Assignments := should be used not =.
2. D cannot be used as a parameter, it is the differentiation operator on functions as in D(sin). Use D1 or similar.
3. List brackets [ ] cannot be used instead of parentheses ( ).

Thus I tried:

restart;
A:=0.1;
D1:=0.19;

 eqns:={y^1.5-9/8*B^0.5*y+3/4*x*y^0.5-3/4*(Pi/3)^0.5*(1+(1+3/4*(Pi/3)*(B*(1+x^2))^0.5)*x^2)/(1+x^2)^0.5=0, Pi/2-3*D1=8/3*y^3+x*y^2+4*B^0.5*(1/3*B*y^1.5-2/45*B^2.5-3/2*y^2.5)-9/2*B*(1/3*B*y-3/40*B^2-3/4*y^2)-3*x*B^0.5*(1/3*B*y^0.5-1/7*B^1.5-1/2*y^1.5),B=A/(1+x^2)};

eq23:=subs(eqns[1],eqns[2..3]); #Thereby eliminating B.

fsolve(eq23,{x,y}); # Returns unevaluated, i.e. no solution found by fsolve.

plots:-implicitplot(eq23,x=-5..5,y=-5..5,gridrefine=1);

It appears that there is no solution for (x,y) in the square ,x=-5..5,y=-5..5.
###################
Changing the parameter D1 to 0.05 e.g. helps. Here is an animation. Notice that D1 is NOT assigned a value.

restart;
with(plots):
A:=0.1;
#D1:=0.19;

 eqns:={y^1.5-9/8*B^0.5*y+3/4*x*y^0.5-3/4*(Pi/3)^0.5*(1+(1+3/4*(Pi/3)*(B*(1+x^2))^0.5)*x^2)/(1+x^2)^0.5=0, Pi/2-3*D1=8/3*y^3+x*y^2+4*B^0.5*(1/3*B*y^1.5-2/45*B^2.5-3/2*y^2.5)-9/2*B*(1/3*B*y-3/40*B^2-3/4*y^2)-3*x*B^0.5*(1/3*B*y^0.5-1/7*B^1.5-1/2*y^1.5),B=A/(1+x^2)};
eq23:=subs(eqns[1],eqns[2..3]);

animate(implicitplot,[eq23,x=0..1,y=0..1.2,gridrefine=1],D1=0.05..0.1);

fsolve(eval(eq23,D1=0.07),{x,y});

######## Point of tangency can be found ##########
implicitdiff~(eq23,y,x): #Works from Maple 13 and up
map(implicitdiff,eq23,y,x): #Works in all versions
eqd:=`=`(op(%)):
fsolve(eq23 union {eqd},{x,y,D1});

I removed ``:

for i from 1 to 3 do
  file_name:=test||i;
  writedata[APPEND](cat("G:/",file_name),A);
end do;

Notice that the 3 files will contain the same data. Was that intended?

In the Programming Guide, Appendix 1, under "Internal Representations of Data Types", you find in the description of SUM the following statement:
"Simple products, such as a*2, are represented as SUM structures."

A little bit earlier it says "Maple syntax: expr1 * factor1 + expr2 * factor2 ...", where factor1, factor2 are numerical.
Thus I think with that description you can think of 3*a as a SUM with one term.
With two terms:
dismantle(3*a+2*b);

@John Starrett

Your problem is of course linear in the coefficients so I tried the straigtforward linear algebra approach:
#The equations to solve:
eqsA:=[seq(eval(Poly1,{x=XA[i],y=YA[i],z=ZA[i]})=WA[i],i=1..33)]:
eqsB:=[seq(eval(Poly1,{x=XB[i],y=YB[i],z=ZB[i]})=WB[i],i=1..33)]:

A,a:=GenerateMatrix(eqsA,[Poly1Coeff]);
ares:=<Poly1Coeff> =~LinearSolve(Transpose(A).A,Transpose(A).a);
#So the result corresponding to FA would be
eval(Poly1,convert(ares,list));

B,b:=GenerateMatrix(eqsB,[Poly1Coeff]);
bres:=<Poly1Coeff> =~LinearSolve(Transpose(B).B,Transpose(B).b);
#So the result corresponding to FB would be
eval(Poly1,convert(bres,list));
#Notice that the determinants are rather large, thus their inverses are rather small:
Determinant(Transpose(B).B);
Determinant(Transpose(A).A);

Therefore I tried raising Digits significantly: to 30. That helped on your FA and FB.

You got a boundary value problem on the interval eta=0..10.
About which point do you want to linearize? And indeed for what purpose?
Writing your equation as a first order system in the usual way, e.g. by using

res:=DEtools[convertsys]([ode1],{},f(eta),eta,Y,YP);

you see that what you need to linearize is the equation for YP[4].

Linearization would have to take place about a point [ Y0[1],Y0[2],Y0[3], Y0[4] ] and could be done with mtaylor:

F:=subs(res[1],YP[4]);
mtaylor(F,[seq(Y[i]=Y0[i],i=1..4)],2);

More conventionally you would use VectorCalculus:-Jacobian on the list of right hand sides:

L:=rhs~(res[1]);
J1:=VectorCalculus:-Jacobian(L,[seq(Y[i],i=1..4)]);
J:=unapply(J1,[seq(Y[i],i=1..4)]): #As a function
J(f0,0,0,0); #At a rather special point: (f0,0,0,0) is an equilibrium for all f0<>0., but unstable.

In your example you have a boundary value problem, thus the parameters option to dsolve is not available. However you could do:

restart;
dsys3 := {H*(diff(f1(x), x, x, x, x))+2*(diff(f1(x), x, x))+2*H*x^2*f1(x) = 10, f1(0) = 0, f1(1) = 0, ((D@@1)(f1))(0) = 0, ((D@@1)(f1))(1) = 0};

#I didn't change all your optional argument, I only left out the range option as it is not needed and I changed array to Array.
#As a trial example I'm evaluating for H=0.12345:

dsol5 := dsolve(eval(dsys3,H=.12345), 'maxmesh' = 500, numeric, abserr = .1, output = Array([.5]));

#If you are going to do this for a lot of H-values you can make a procedure:

p:=proc(H1) dsolve(eval(dsys3,H=H1), 'maxmesh' = 500, numeric, abserr = .1, output = array([.5])) end proc;
#Example of the use:
p(.8);

If you actually want f1(0.5) as a function of H you can do:

q:=proc(H1) local dsol,R;
   dsol:=dsolve(eval(dsys3,H=H1), 'maxmesh' = 500, numeric, abserr = .1, output = Array([.5]));
   R:=dsol[2,1];
   R[1,2]
end proc;
q(.456);
plot(q,0.2..1); #You may want to change abserr to something smaller

There is no reason for all the splitting:

restart;
epsilon:=5:Delta1:=4:Delta2:=4:  
dsys :={diff(x(t),t)=-I*Delta1*x(t)+y(t)+epsilon, diff(y(t),t)=-I*Delta2*y(t)+x(t)*z(t), diff(z(t),t)=-2*(conjugate(x(t))*y(t)+conjugate(y(t))*x(t))};

res:=dsolve(dsys union {x(0)=0,y(0)=0,z(0)=-1/2},numeric);

plots:-odeplot(res,[[t,Re(x(t))],[t,z(t)]],0..2);

To get you started:

restart;
sys:=diff(x(t),t,t)=a*x(t)+b*y(t),diff(y(t),t,t)=c*x(t)+d*y(t);
ics:=x(0)=1,D(x)(0)=0,y(0)=0,D(y)(0)=1;
#Solving in general:
res:=dsolve({sys,ics});
#Solving with the particular choice of parameters:
res2:=dsolve(eval({sys,ics},{a=1,b=1,c=1,d=1}));
# Extracting the x-solution from res2:
X2:=subs(res2,x(t));
# Evaluating the general result for the particular choice of parameters:
eval(res,{a=1,b=1,c=1,d=1}); #Division by zero: special case
#Extracting x:
X:=subs(res,x(t)):
#Trying limit as d-> 1:
X1:=limit(eval(X,{a=1,b=1,c=1}),d=1);
X1:=combine(expand(X1)) assuming t::real;
# Is X1=X2?
X1-X2;
simplify(%);
#Yes.
Conclusion: In res you have the generic case, i.e. the solution except for special values of a,b,c,d (actually the ones for which the determinant of the matrix Matrix([[a,b],[c,d]]) is zero).
Taking limits recovers the special case. Better to start with the special case, if that is all you want to look at.
##############
# I notice that you have the tag Laplace Transform.
So you can try:
resL:=dsolve({sys,ics},{x(t),y(t)},method=laplace);
allvalues(resL);
eval(%,{a=1,b=1,c=1,d=1});
simplify(%); #The division by zero hasn't gone away, but has to be dealt with similar to the above.

Finally, I think I really know what is causing problems for all output options except output=Array and output=piecewise.

If you plot before saving, as in

plots:-odeplot(S1,[x,z(x)],0..1);

then the version read in a new session doesn't work.