@Rouben Rostamian  That is very nice indeed.
A comment:

1. It could be mentoned that while x(t) = xtr+xss is a solution to de neither xtr or xss is a solution to the de.
2. I do believe,however, that there must be a solution defined on the whole real axis. That solution is periodic. All other solutions will tend to either -infinity or +infinity as t -> -infinity.

## The argument for part of point 2 is that the the horizontal strip { (t,x) |  10/11 < x < 10/9 } is positively invariant. In fact the right hand side of the ode f(t,x) is negative for any x > 10/9 and positive for x < 10/11.
Outside the strip the solutions will grow exponentially backwards in time. The ones above the strip will go to +infinity, those below to - infinity.
Consider the x-interval  I = (10/11, 10/9). Starting at any arbitrary time t0 with x(t0) in that interval there are 3 possibilities when going backwards in time: Either the solution leaves the strip upwards to infinity (case 1), downwards to minus infinity (case 2), or stays in the strip for all times (case 3).

Using existence and uniqueness of the initial value problem for this ode we see that the solutions belonging to cases 1 and 2 correspond to non-overlapping open subintervals of the interval I.
The interval of x-values corresponding to case 3 must be closed and not empty, possibly consisting of just one point.
This still leaves open the possibilty that this last interval consisting of solutions that stay in the strip for all negative as well as positive times consists of more than one point.

In our present example surely it seems that there is only one point.

Notice that that x-point will depend on the chosen value of t0. When t0 = 0 is chosen it is something like x0 = 1.00997547504228. Because of the obvious instability in the negative direction solving numerically from (0, x0) will not give us the solution in that direction:

How do you define an equibrium point for a non-autonomous ode?

If you mean a constant solution, i.e. x(t) = k, where k is a constant, then it is trivial to see that your ode doesn't have such a solution.
By inserting x(t) = k into x'(t)=-(1+0.1*sin(10*t))*x(t)+1  you get 0 = -(1+0.1*sin(10*t))*k+1 = -k - 0.1*k*sin(0.1*t) +1

Since this has to hold on an interval t = a..b k must be zero. Thus we get 0 = 1, a contradiction!

What you may be able to show is that all(?) solutions tend towards some common periodic solution x0(t) as t -> infinity, i.e.
x(t)-x0(t) -> 0 as t -> infinity.

It is easy to give convincing numerical evidence of that claim.

@AHSAN I would recommend Maple's own manuals. There is a User Manual and a Programming Guide.
Go to https://www.maplesoft.com/documentation_center/

They have manuals from previous versions too.

Furthermore, I would strongly recommend that you use Maple Input (1D math input) instead of 2D and also Worksheet Mode.

To make that change in Maple go to the menu item Tools. Then Options/Display/Input display. Here choose Maple Notation.

After that and still with the Options window open go to Interface/Default format for new worksheets. Here choose Worksheet.

Now click on Apply Globally. Don't worry, the changes can easily be undone.
The changes only affects new worksheets. So try opening a new worksheet and start writing some Maple code on an input line marked with >.
###########################################

I have added a link to a worksheet in the format recommended here.

The worksheet includes the code shown above, but also conversion of the second order ode to a first order system in the usual way. Furthermore this system is solve exactly as well as numerically by a shooting method.

MaplePrimes21-08-05_ode_bvp.mw

@AHSAN You ignored my correction of your syntax error. Please see above.
Since you refer to a paper, where they use numerical solution, and not a textbook teaching you about bvps, there seems to be absolutely no point in not using the exact method given above.
I was beginning to think that your instructor had given you this problem as an exercise and wanted it done in the way you describe.

That not being so just use the result sol1 or sol (sol = sol1) as found also by Maple 2019.2.

That result is:

u(y) = -1/2*y/beta + 1/12*(4*beta*p*y + 1)^(3/2)/(beta^2*p) + 1 + 1/2*h/beta - 1/12*(4*beta*h*p + 1)^(3/2)/(beta^2*p);

From that you can compute all the values you want.
Incidentally, why didn't you update your Maple 2019.0 to 2019.2?

@vs140580 I don't have Excel on my computer, so I cannot help you with that.
I'm assuming that you know of the ExcelTools package. It has Export and Import to and from Excel.

@nm Using the workaround (as it seems to be) by Joe Riel, i.e. using identical(x)^2 instead of identical(x^2) and using *& with select as in my answer, you also have to be aware of the non-commutativity of *&.

restart;
expr:=x^2*A+B;
select(type,expr,identical(x)^2 &* anything);
select(type,expr,anything &* identical(x)^2);
select(type,expr,{anything &* identical(x)^2,identical(x)^2 &* anything}); # safer

With indexed names:

restart;
expr:=x^2*A[2]+A[1];

select(type,expr,identical(x)^2 &* anything);
select(type,expr,anything &* identical(x)^2);
select(type,expr,{anything &* identical(x)^2,identical(x)^2 &* anything});

@nm You are right!

@Sradharam Do you mean to say that the boundary conditions for theta are theta(0)=1 and theta(infinity)=0 instead of what you actually wrote?

If your answer is YES then in my code just make the following changes:
1. In the procedure Q change TH(inf) -1 to TH(inf).
2. Use this th1 interval th1=-0.1325..-0.1275 in the animation and the same interval in the 3 solving methods.

Make no other changes.
You get:

params := {f1 = -0.9272464891, th1 = -0.1295294672}

and the graphs:

@Sradharam I admit that the code is in need of some explanation.
Here I will just point to a few essential features.

1. Using the parameters option in dsolve ( here parameters=[f1,th1] ) sets up the result (here named respar) as a complete recipe for computation before any values are actually given to those parameters. This avoids having to call dsolve for every new pair of actual values for the parameters. To insert values for the parameters you just do
respar( parameters = [xx, yy]) or more explicitly respar( parameters = [f1=xx, th1=yy]), where xx and yy are actual numbers.

2. Using output = listprocedure instead of the default output = procedurelist makes it easy to extract the procedures you will actually use. In our case they are called F and TH. They compute values for f(xi) and theta(xi), respectively, when given a concrete value for xi. We need F(inf) = 0 and TH(inf) -1 = 0, where inf is here taken to be 10 later.

3. The procedure Q produces a list of the values F(inf) and TH(inf) -1 on input of concrete values for f1 and th1.
To make it possible to solve for f1 and th1 using either of the 3 methods I used I have made the procedures q1 and q2 return the single values F(inf) and TH(inf) -1, respectively. The option remember in Q implies that a call to q1 will make it very cheap to answer the call to q2. Both will be taken from the same call to Q.

4. Since respar( parameters = [xx, yy]) sets the parameters to the values xx and yy the parameters are reset many times when either of the 3 solution methods are used. Since we don't necessarily know what values are set last, the cautious user will set the parameters returned by either of the 3 methods if respar (or F or TH) is going to be used with those parameters. I used respar subsequently to plot using odeplot. I set respar at the result from fsolve.

5. I'm using the default rkf45 dsolve method as mentioned. This is a variable step method. It adjusts the step taken in order to make sure that the absolute and relative error bounds are satisfied. The classical RK4 method is a fixed step method, i.e. you choose the stepsize in advance. There is no error control.
Another drawback for RK4 is that the parameters option is not available for that method.

@Sradharam As Tom Leslie points out, you need to give us the value of Pr.

I have taken it to be 0.01, and I use alpha=2, k=1.5.
I'm using the same shooting method I used above, but it is now extended to two variables.
The numerical method used by dsolve is the default rkf45 (rkf=Runge-Kutta-Fehlberg).
 

restart;
odeSys := f(xi)+3*diff(g(xi), xi)-xi*diff(f(xi),xi)=0,
            f(xi)^2+3*g(xi)*diff(f(xi),xi)-xi*f(xi)*diff(f(xi),xi)=3*diff(f(xi),xi,xi)+18*k*diff(f(xi),xi)^2*diff(f(xi),xi,xi);
  bcs := f(0)=alpha, f(infinity)=0, g(0)=0;
ode3 := theta(xi)*f(xi) + g(xi)*diff(theta(xi), xi) - 1/3*xi*f(xi)*diff(theta(xi), xi) = diff(theta(xi), xi, xi)/Pr;
bcs_theta:=theta(0) = 1, theta(infinity) = 1;
####
respar:=dsolve(eval({odeSys,ode3,f(0)=alpha,g(0)=0,D(f)(0)=f1,theta(0)=1,D(theta)(0)=th1},{alpha=2,k=1.5,Pr=.01}),numeric,
                parameters=[f1,th1],output=listprocedure,abserr=1e-12,relerr=1e-10);
F,TH:=op(subs(respar,[f,theta](xi)));
Q:=proc(f1,th1,inf) option remember;
#if not [f1,th1,inf]::list(realcons) then return 'procname(_passed)' end if;
   respar(parameters=[f1,th1]);
   [F(inf),TH(inf)-1]
end proc;
q1:=proc(f1,th1,inf) if not [f1,th1,inf]::list(realcons) then return 'procname(_passed)' end if;
       Q(_passed)[1] 
end proc;
q2:=proc(f1,th1,inf) if not [f1,th1,inf]::list(realcons) then return 'procname(_passed)' end if;
       Q(_passed)[2] 
end proc;
###
### Some exploratory plotting:
plot([q1(f1,-0.025,10),q2(f1,-0.025,10)],f1=-2..0);
plots:-animate(plot,[[q1(f1,th1,10),q2(f1,th1,10)],f1=-1..-0.9,-0.3..0.3],th1=-0.03..0);
### From the animation we see that the values for theta must lie in the range:
th1=-0.025..-0.02375; # sol1
### in order to find the relevant solution (there are 3 as before).
### I find the values for f1 and th1 in 3 ways although one is enough.
### Skip the first if you don't have DirectSearch yet.
with(DirectSearch);
SolveEquations([q1(f1,th1,10),q2(f1,th1,10)],{f1=-0.95..-0.92,th1=-0.025..-0.02375});
Optimization:-LSSolve([q1(f1,th1,10),q2(f1,th1,10)],f1=-0.95..-0.92,th1=-.025..-0.02375);
params:=fsolve([q1(f1,th1,10),q2(f1,th1,10)],{f1=-0.95..-0.92,th1=-.025..-0.02375});
respar(parameters=convert(params,list)); # Making sure parameters are set at params.
plots:-odeplot(respar,[[xi,f(xi)],[xi,g(xi)],[xi,theta(xi)]],0..10,thickness=3,legend=[f,g,theta]);
plots:-odeplot(respar,[xi,theta(xi)],0..10);

Here are the graphs:

Just downloaded your worksheet and executed the whole, but got (in Maple 2021.1):

Error, (in plottools:-rotate) invalid plot structure: wheel

What is p1 in the code for wheel in the line
 

p[i++] := translate(p1, 0,0,-h):

?

Added: It appears that d1 should be %.

Also in the code for spokes wheel should be wheel().

@Carl Love In Maple 2021.1 numerical integration of KummerM doesn't work with method=_d01ajc. Thus h[1] [1] doesn't get evaluated.in the double loop. I have submitted an SCR.
 

restart;
K:=KummerM(-0.0847656450, 5., 10.33906258*r^2);
int(K,r=0..1,numeric, method = _d01ajc);

@mmcdara To your first question: That is what I thought of first. You are right, your shorter version works too.

The type And(relation,satisfies(r->lhs(r)=t)) is satisfied by expressions that are both of type relation and also of the special type 'satisfies'. The argument to 'satisfies' is a procedure that returns true (in this example) if the relation (here t < something, or x(t)<something) has left hand side equal to t.
See type/satisfies.

@itsme You may be right. The acceptance of vector input is fairly new.

By the way there is another inconsistency: dsolve with method=laplace or type=series will not work without the variable(s) as a second argument. That I have been told is due to different people writing those parts.
 

restart;

sys:=Vector([diff(x(t),t),diff(y(t),t)]) = Matrix([[1,2],[0,3]]).Vector([x(t),y(t)]);
init:=<x(0),y(0)> = <x0,y0>; # vector equation
### I think the form of the following outputs are quite fine:
dsolve(sys);         # output set of eqs
dsolve({sys,init}); # output set of eqs
## Giving the second argument as a vector forces vector output:
dsolve({sys,init},<x(t),y(t)>); # output vector equation
dsolve(sys,<x(t),y(t)>);         # output vector equation

### The next outputs are also fine.
### The second argument, however, MUST be present when using method = laplace or (type=)series

dsolve(sys,[x(t),y(t)],method=laplace);   # output set of eqs
dsolve(sys,[x(t),y(t)],series);                   # output set of eqs
dsolve(sys,<x(t),y(t)>,method=laplace); # output vector equation
dsolve(sys,<x(t),y(t)>,series);                 # output vector equation

@Carl Love In the help page for seq the following two illustrations are given. Here adapted to the present discussion.
 

restart;
S:=NULL;
old:=x[4];
for x[4] from 1 to 10 do S:=S,u(x[4]) end do;
x[4]:=old;
S;
restart;
S:=NULL;
old:=x[4];
for x[4] in [1,2,3,4,5,6,7,8,9,10] do S:=S,u(x[4]) end do;
x[4]:=old;
S;

Here we get the assignment 11 to x[4] in the first and 10 to x[4] in the second.