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

Change the dsolve line to:

dsol[i] := CodeTools:-Usage(dsolve(dsys, numeric, continuation = lambda));

You mentioned method RK45 (rkf45?) in your title. Your problem, however, is a boundary value problem and is handled by `dsolve/numeric/bvp`. Method rkf45 is used for initial value problems only.

You can solve the problem by the textbook method of Fourier series.
The first step is a well known trick which is used to reduce the problem to a homogeneous one.

restart;
pde:=diff(u(x, t), t, t)-c^2*(diff(u(x, t), x, x))-cos(x) = 0;
IBC :={u(0, t) = A, u(1, t) = B, u(x, 0) = 0, (D[2](u))(x, 0) = sin(x)};
eval(pde,u(x,t)=a*x+b+cos(x)/c^2);
sol1:=u(x,t)=a*x+b+cos(x)/c^2;
eval(rhs(sol1),x=0)=A;
eval(rhs(sol1),x=1)=B;
solve({%%,%},{a,b});
sol2:=eval(sol1,%);
pde2:=simplify(eval(pde,u(x,t)=v(x,t)+rhs(sol2)));
pdsolve(pde2,HINT=X(x)*T(t));
dsolve(diff(X(x), x, x) = -k^2*X(x));
solve(sin(k)=0,k,allsolutions);
#k = n*Pi
Tsol:=dsolve(diff(T(t),t,t) = -c^2*n^2*Pi^2*T(t));
S:=Sum(sin(n*Pi*x)*(a(n)*cos(Pi*n*c*t)+b(n)*sin(Pi*n*c*t)),n=1..infinity); #Solution for v(x,t)
eval(S,t=0)=-rhs(sol2); #u(x,0)=0
res_a:=a(n)=2*int(rhs(%)*sin(n*Pi*x),x=0..1) assuming n::posint;
diff(S,t); #We can ignore the time independent sol1
eval(%,t=0)=sin(x);
res_b:=b(n)=1/Pi/n/c*2*int(sin(x)*sin(n*Pi*x),x=0..1) assuming n::posint;
S1:=eval(S,{res_a,res_b});
res:=u(x,t)=rhs(sol2)+S1;
plots:-animate(plot,[eval(rhs(res),{A=5,B=9,c=1,infinity=50}),x=0..1],t=0..10,frames=100);

As I said in my comment, you need to rewrite the boundary conditions to get rid of f'''(0) or in the first order systems version v'(0).
Here I do that for the first order system. I do not use a shooting method.
#####
restart;
varsigma := 10;
a := 1;
b := -1;
ss := 2;
Pr := 1;
lambda := 0;
sigma := -2;
ODE := {(diff(t(eta), eta))/Pr+f(eta)*t(eta)-u(eta)*s(eta) = 0, diff(v(eta), eta)+f(eta)*v(eta)-u(eta)^2+lambda*s(eta) = 0, diff(f(eta), eta) = u(eta), diff(s(eta), eta) = t(eta), diff(u(eta), eta) = v(eta)};
## Your boundary conditions:
bcs:={u(varsigma) = 0, s(varsigma) = 0,f(0) = ss, u(0) = sigma+a*v(0)+b*(D(v))(0), s(0) = 1};
## We need to get rid of D(v)(0).
## So we use the ode for v:
ode_v:=diff(v(eta), eta)+f(eta)*v(eta)-u(eta)^2+lambda*s(eta) = 0;
eval(convert(ode_v,D),eta=0);
dv0:=solve(%,{D(v)(0)});
## The rewritten boundary conditions:
bcs2:=eval(bcs,dv0);
res:=dsolve(ODE union bcs2,numeric);
plots:-odeplot(res,[[eta,f(eta)],[eta,s(eta)]],0..varsigma);

The first order systems approach gave answers right away. The original version you posted gave some problems, I guess because of the one nonlinear boundary condition. This condition is also present in the first order version, but just happened not to cause a problem.
I tried using the result from the first order system as an approximate solution for the original. It was comforting that that worked fine:
ode1:=diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta,eta)-diff(f(eta),eta)^2+lambda*theta(eta)=0;
ode2:=1/Pr*diff(theta(eta),eta,eta)+f(eta)*diff(theta(eta),eta)-diff(f(eta),eta)*theta(eta)=0;
## We need to rewrite this:
bcsE:={f(0)=ss,D(f)(0)=sigma+a*(D@D)(f)(0)+b*(D@@3)(f)(0),theta(0)=1,D(f)(varsigma)=0,theta(varsigma)=0};
eval(convert(ode1,D),eta=0);
f30:=solve(%,{(D@@3)(f)(0)});
bcsE2:=eval(bcsE,f30); #The rewritten boundary conditions
## Now we use the result from the first order system:
resM:=dsolve(ODE union bcs2,numeric,output=mesh); #output=mesh is easiest to deal with right now.
## We need to change names and such:
resM2:=subs(s(eta)=theta(eta),t(eta)=diff(theta(eta),eta),u(eta)=diff(f(eta),eta),v(eta)=diff(f(eta),eta,eta),resM);
## Now using resM2 as an approximate solution:
sol:=dsolve({ode1,ode2} union bcsE2,numeric,approxsoln=resM2);
plots:-odeplot(sol,[[eta,f(eta)],[eta,theta(eta)]],0..varsigma); # Same as shown above
##Try first without an approximate solution. You will see problems, which I'm sure is because of the one nonlinear boundary condition.
##############
## Exploring the effect of rewriting the nonlinear boundary condition in the latter formulation:
eqE:=select(has,bcsE2,D(f)(0));
#Choosing to isolate f '' (0):
d20a:=eval(solve(eqE,{(D@D)(f)(0)}),f(0)=ss); # Also explicitly putting f(0)=ss.
d20b:=solve(eqE,{(D@D)(f)(0)}); # Leaving to dsolve to use f(0)=ss.
bcsE3a:=remove(has,bcsE2,D(f)(0)) union d20a;
bcsE3b:=remove(has,bcsE2,D(f)(0)) union d20b;
## Now solving the system with these two formulations of the boundary conditions. We use the same approximate solution in the two cases.
## Using bcsE3a:
resE3a:=dsolve({ode1,ode2} union bcsE3a,numeric,approxsoln=[f(eta)=ss*exp(-2*eta)+ss,theta(eta)=exp(-2*eta)],initmesh=2048);
plots:-odeplot(resE3a,[[eta,f(eta)],[eta,theta(eta)]],0..varsigma); pa:=%:
## Using bcsE3b:
resE3b:=dsolve({ode1,ode2} union bcsE3b,numeric,approxsoln=[f(eta)=ss*exp(-2*eta)+ss,theta(eta)=exp(-2*eta)],initmesh=2048);
plots:-odeplot(resE3b,[[eta,f(eta)],[eta,theta(eta)]],0..varsigma); pb:=%:
plots:-display(pa,pb);  #Surprise! Two solutions.

#Finally (I hope) I tried my own shooting procedure and was without problem able to reproduce the two different solutions.

You don't tell us what f__r and f__h are, but I made up a simple example, where I replaced your large number 9999999 by 'undefined':

#First using 9999999:
restart;
q:= proc (x,y,z)
  if x>0 then x^2+y^2+z^2-1 else 9999999 end if
end proc;
plots:-implicitplot3d('q(x,y,z)',x=-2..2,y=-2..2,z=-2..2);
#Now using 'undefined':
restart;
q:= proc (x,y,z)
  if x>0 then x^2+y^2+z^2-1 else undefined end if
end proc;
plots:-implicitplot3d('q(x,y,z)',x=-2..2,y=-2..2,z=-2..2);