You are clearly using Douglas Meade's shoot procedure.
You should not use period (.) for ordinary multiplication. Use *.
You have not solved for the derivatives.
To correct these two problems do:
ODE1:=eval(ODE,`.`=`*`);
dfs:=indets(ODE1,specfunc(diff));
ODE2:=solve(ODE1,dfs);
S := shoot(ODE2, IC, BC, FNS, [alpha1 = 1.425, alpha2 = .425, beta1 = -1.31, beta2 = 1.00, beta3 = 1.29]);
odeplot(S, [eta, fp(eta)], 0 .. 1.84,view=-100..100);
#You will observe a singularity problem close to eta=1.84.
odeplot(S, [eta, g(eta)], 0 .. 1.84,view=-10..1);

Maybe like this:

f:=cos(n*t)+n; #Example
S:=seq(plots:-spacecurve([t,n,f],t=0..2*Pi,labels=[t,'n','f']),n=0..10):
plots:-display(S,orientation=[-10,75,5]);

One of several undocumented "methods"  _RETURNVERBOSE reveals the problem.
As you noticed, the problem is in the first integration (over x). The answer is zero if we simplify to y>-1,y<1,y<>0, to avoid unnecessary complication:
int(Heaviside(-x^2-y^2+1), x = -1 .. 1) assuming y<1,y>-1,y<>0; #Answer 0
## Now try
int(Heaviside(-x^2-y^2+1), x = -1 .. 1,method=_RETURNVERBOSE)  assuming y<1,y>-1,y<>0;
## You will see that piecewise, ftoc, ftocms, and distribution report successes. The remaining fail.
distribution reports 0, the other 3 report 2*sqrt(-y^2+1)
## The distribution result is obviously wrong. So try any of the other 3:
int(Heaviside(-x^2-y^2+1), x = -1 .. 1,method=ftoc) assuming y<1,y>-1,y<>0;

#The error is also present in Maple 2016.
There are other "methods" besides _RETURNVERBOSE, e.g. the following
int(Heaviside(-x^2-y^2+1), x = -1 .. 1,method=_RETURNALL)  assuming y<1,y>-1,y<>0;

## The double integral then could be done like this:
int(int(Heaviside(-x^2-y^2+1), x = -1 .. 1,method=ftoc),y=-1..1);

#NOTE.  I have submitted an SCR suggesting that some of the undocumented integration "methods" mentioned above be documented.

I'm using Maple 2016 here (but same results in Maple 12):

u:=sqrt(sqrt(9)*sqrt((1+(b+1)^2*c^2+((10/3)*b-2)*c)*(1+(b+1)^2*c^2+(2*b-2)*c))+3+(3*b^2+6*b+3)*c^2+(8*b-6)*c);
## Because sqrt in Maple gives the principal square root u > 0 iff u^2 > 0:
u2:=u^2;
coulditbe(u2>0); #true in Maple 2016 and in Maple 12
coulditbe(u>0); #true in Maple 2016 and in Maple 12
plot3d(u2,b=0..5,c=0..5);#[0,5] Obviously there are plenty of values of (b,c) in the square [0,5]x[0,5] for which u2 (thus also u) is positive.
#####################################
In fact your expression is positive for all real nonzero b and c for which the square root in u2 is real (thus >=0).
With:
v,h:=selectremove(hastype,u2,radical);
# u2 = 0 implies
v^2=h^2;
expand((lhs-rhs)(%)); # Answer -4*b^2*c^2 which is only zero if b=0 or c=0.
#By plotting as above we can see that there are points in all 4 quadrants for which u2 > 0.
#Thus u2 > 0 for all (b,c) with b<>0 and c<>0 since u2 is continuous in (b,c).
#What remains is to ask if the quantity under the square root in u2 given by
w:=(1+(b+1)^2*c^2+((10/3)*b-2)*c)*(1+(b+1)^2*c^2+(2*b-2)*c);
# is nonnegative for all real (b,c).
# That appears to be the case:
plot3d(w,b=-5..5,c=-5..5);
# and minimize confirms it:
minimize(w);

Try
res:=dsolve(ode);
#The output is a sequence of two results on implicit form and both involving an integral, which cannot be computed #unless n is specied. Even with n=0 the result doesn't look simple:
value(eval([res],n=0));

This produces what you want:

series(cm^(-1),cm);

There may be other less strange ways of doing this.

First of all you can certainly do this numerically, but you can also find the solution using fourier series.

restart;
pde:=  diff(u(x,t),t)=diff(u(x,t),x,x);
## First numerically:
res:=pdsolve(pde,{u(0,t)=t,u(1,t)=5*t,u(x,0)=x},numeric);
res:-animate(t=10,frames=50);
###
Now using fourier series and a trick, which reduces the problem to a problem with homogeneous boundary conditions.
## The trick is this. Consider this function:
V:=u(x,t)=x*5*t+(1-x)*t+a*x^3+b*x^2+c*x;
# where a,b,c are to be determined. We want V to be a solution to pde and satisfy the boundary conditions:
eq1:=eval((lhs-rhs)(pde),V);
res1:=solve(identity(eq1,x),{a,b}); #Requirement for pde to be satisfied.
eval(V,x=0); #Satisfies boundary condition at x=0.
eval(V,{x=1} union res1);
res2:={c=-7/6}; #Requirement for boundary condition at x=1.
sol1:=eval(V,res1 union res2); # The solution for V.
##Just to be sure we check the answer:
eval(sol1,x=0);
eval(sol1,x=1);
eval((lhs-rhs)(pde),sol1);
## So OK
##Now we write the solution to the whole problem as u(x,t)=v(x,t) + u0(x,t), where v(x,t) is sol1 and u0(x,t) is to be found by using fourier series.
## For u0(x,t) we must use the following initial value:
ic:=u(x,0)=x-eval(rhs(sol1),t=0);
# Now Maple 2016 can find u0(x,t):
solH:=pdsolve({pde,u(0,t)=0,u(1,t)=0,ic});
## For convenience we get rid of _Zn~ (n some integer) (this is not necessary).
z:=op(indets(solH,name) minus {Pi,x,t,infinity});
sol:=subs(z=n,rhs(solH)+rhs(sol1));
plots:-animate(plot,[eval(sol,{infinity=100,Sum=add}),x=0..1],t=0..10,frames=50);

It is comforting to observe that the numerically found animation (not shown here) agrees well with the fourier version.

Try

f(31536000); # No period, i.e. input not a float
## Then
evalf(%);
                 2.718281785

Well, your equation is most likely constructed by your instructor (wild guess).

The answer is that y(x) = exp(x) is a solution. (Are there other solutions?)

Check using the reformulation by vv:

ide:=x^2*(diff(y(x), x, x))+50*x*(diff(y(x), x))-35*y(x) - (1-exp(x+1))/(x+1)-(x^2+50*x-35)*exp(x)-int(exp(x*t)*y(t),t=0..1);

simplify(eval(ide,y=exp)); #result 0.

Now at least you have something with which to compare your numerically found solution.

There is no chance of getting an exact answer. But you can get a numerical approximation to any of the infinitely many roots by using fsolve:

Digits := 20;
L := 20*R;
eq:=tan(sigma*L/(2*R))+tanh(sigma*L/(2*R))=0; #No names except sigma
plot(rhs(eq),sigma=-1..1,discont=true);
fsolve(eq, sigma = 0.5); #Finding the root near 0.5

Why don't you want to use Maple's numeric bvp solver? It uses a finite difference method.

I gave a solution in this link:

http://mapleprimes.com/questions/212370-About--Nonlinear-ODE-System-Again

For convenience I copy it here:

You can do as follows.
The idea is to start with the approximate solution provided by you. But in the next step the previously found solution is used for the new value of s1. For convenience the results are kept in vectors.

Sigma:=Vector(datatype=float);
S1:=Vector(datatype=float);
AP:=[omega2 = 0.5e-1, s(x) = (56330/77299)*x^6-(25403/9094)*x^5+(1214873/309196)*x^4-(182218/77299)*x^3+(1/2)*x^2, g(x) = -(32119348934425/2772462826064)*x^5+(133402088597045/5544925652128)*x^4-(35403636020221/2772462826064)*x^3-(1/2)*x^2+x];
i:=0:
for s1 from 0.001 to 0.06 by 0.001 do
   try
     res1:=dsolve(dsys4 union {s(1) = s1}, numeric, initmesh = 1024, approxsoln = AP, abserr = 0.1e-4,maxmesh=8192);
     AP:=res1;
     i:=i+1;
     Sigma(i):=eval(sqrt(omega2(1))/(L*sqrt(m/E_a)),res1(1));
     S1(i):=s1;
   catch:
     print(s1);
     print(lasterror);
   end try;
end do; 
plot(Sigma,S1,labels=[sigma,s(1)]);

You have reformulated your earlier question
http://mapleprimes.com/questions/213849-Plot-A-Sequence--Of-Functions-Between-2-Colors

Here is an example that resembles the present situation:

n:=10:
T:=[$1..n+1];
S:=[sin(i*x)$i=1..n];
P:=seq(plot(S[i],x=T[i]..T[i+1],color=`if`(i::odd,red,black)),i=1..n):
#plots:-display(P,size=[1800,default]); # remove the size option for Maple 14. It only works in recent versions of Maple.
plots:-display(P); #For older versions

Here is an example which should work whether n is even or odd:

n:=11:
L:=[seq(x^i,i=1..n)];
plot(L,x=-1..1,color=op~([[red,black]$ceil(n/2)]));

I don't think that fsolve uses an interval unless such is provided. But it does need a starting guess just as Newton's method does. The usual Newton method doesn't require an interval.
If a starting guess is not provided fsolve will find one itself.

To explore try the following where p(x) is just sin(x) (as long as x has type numeric). But as a side effect p prints input and output.

restart;
p:=proc(x) local res; res:=sin(x); printf("input %f output %f\n",x,res); res  end proc;
fsolve(p); #No guess provided
fsolve(p,1000); #Initial guess 1000
fsolve(p, 5); # Initial guess 5
fsolve(p,Pi/2); # Initial guess Pi/2 where the derivative of sin (i.e. cos) is zero
fsolve(p,Pi/2..4); #An interval is provided