@janhardo At the end I wrote:
" After that you could select the whole text line and convert to plain text.
In fact this is all you need to do."
You just select the text line and right click. Choose convert to Plain Text from the the menu coming up.

Then that error disappears, if the worksheet is run again.

It appears to me that the term w*S(t) in the ode for S(t) should be w*C(t) when I compare to the pdf file.
Thus in the notation used by Tom Leslie your system should be:

odeSys:= diff(S(t),t) = z*(1 - P[1] - P[2]) + w*C(t) - x*Q[1]*C(t)*S(t)/N - x*Q[2]*Iota(t)*S(t)/N - v*S(t),
           diff(C(t),t) = z*P[1] + x*Q[1]*C(t)*S(t)/N + x*Q[2]*Iota(t)*S(t)/N - w*C(t) - k*o*C(t) - v*C(t),
           diff(Iota(t),t) = C(t)*k*o - Iota(t)*v + z*P[2];

in which also the correction pointed out by Carl Love has been made.
Then we get that the increase in the total:  Total(t) = S(t) + C(t) + I(t) satisfies a very reasonable ode:
 

collect(add(odeSys[j],j=1..3),[z,v]);
odeTotal:=eval(%,Iota(t)=Total(t)-S(t)-C(t));
dsolve({odeTotal,Total(0)=N0});

We find odeTotal := diff(Total(t), t) = z - Total(t)*v, which is also what we get from the odes in the pdf file when z = pi and v = delta.

Notice that Total(t) is not constant in general, so the quantity N (treated as a constant) must be just some convenient number.

@Reshu Gupta For R = 10 you can use e.g.  maxmesh=1024,initmesh=512.
So make use of the options given in the help.

@Reshu Gupta In what sense is the result from dsolve/numeric/bvp inappropriate or not appropriate enough?
If you mean accuracy, you could try raising abserr.

@Reshu Gupta Your given boundary conditions IC (misleading if I = initial and C = condition) are:
IC:={D(D(f))(0)=0, D(f)(1)=0,f(0)=0,f(1)=1,g(0)=0,  g(1)=0}:
What do you expect Maple to do with those? They are given by you as are also the differential equations EQ.
When you say:
" I need to solve it by Runge-Kutta fourth order Method. Kindly help me in the coding of the same. "
Why do you think you need to do that?

@Axel Vogt Your way is nice and simple!
 

rationalize(expr); 
res:=int(rationalize(expr),p); 
simplify(diff(res,p)-expr); # 0

This even works in Maple 8. In fact in Maple 8 rationalize is not necessary.
Maple 8 comes up with the following without assumptions just from res:=int(expr,p);
 

res := 1/2*(2*arctan(p)*a^2+2*arctan(p)*a^2*p^2+4*arctanh(p/(a^2-1)^(1/2))*(a^2-1)^(1/2)*p^2-((p^2+1)^2*(a^2-1))^(1/2)*ln(p^2+1)+4*arctanh(p/(a^2-1)^(1/2))*(a^2-1)^(1/2)+2*((p^2+1)^2*(a^2-1))^(1/2)*ln((p^3*(a^2-1)^(1/2)+p*(a^2-1)^(1/2)-((p^2+1)^2*(a^2-1))^(1/2)*((a-1)*(a+1))^(1/2))/(p^2+1)/(a^2-1)^(1/2))-2*arctan(p)-2*arctan(p)*p^2)/((p^2+1)^2*(a^2-1))^(1/2);
###
simplify(diff(res,p)-expr); # 0

Obviously you didn't give us the full code.

Could you do that or simply upload a worksheet with the whole shebang?

Runge-Kutta methods are used for initial value problems, but you have a boundary value problem.

@Carl Love Exploring this just one step further indicates a bug in Maple 2020.1, Windows 10, June 10 2020 Build ID 1474787.

restart;
u:=1/3*ln(_C1^2 - 2*_C1*x + x^2 + 1) + 1/2*(-2*x + 2*_C1)*arctan(-x + _C1);
simplify(u);

I get
1/3*ln(_C1^2 - 2*_C1*x + x^2 + 1) + 1/3*(-3*x + 3*_C1)*arctan(-x + _C1)

which is correct, but couldn't possibly be intended.
######################

A simpler example:
 

restart;
u:=1/3*g(a) + (x+y)*f(b);
simplify(u);

Result:
1/3*(3*x + 3*y)*f(b) + 1/3*g(a)

PS. This goes all the way back to Maple 2016.2. It is not present in Maple 2015.2.

@Carl Love I need to use simplify directly on that term:
 

restart;
interface(version); # `Standard Worksheet Interface, Maple 2020.1, Windows 10, June 10 2020 Build ID 1474787`
my_sol:=y(x) = arctan(x - _C1)*x - arctan(x - _C1)*_C1 - ln(1 + (x - _C1)^2)/2;
simplify(my_sol); # 2 is still there
map(simplify,rhs(%)); # 2 is gone

@nm dsolve/numeric/bvp comes up with the initial profile n(x) = 0, which obviously is no good since n(x) is a factor in the denominator of the expression for n''(x).

So if 'Vortex' has an idea about the shape of a solution he could try giving an approximate solution in the form
approxsoln=[n(x) = f(x)]
where f(x) is something having that shape.

@AHSAN Since you don't say what it is that needs more explanation I shall make a few comments about what I think may require some elaboration. I will restrict myself to the last and fastest version.

1. The use of unapply. In your own worksheet you have the lines:
 

lambda1 := -3*(7*k^3*sigma^3 + 32*Q*k^2*sigma^2 - 11*k^2*sigma^3 + 54*Q^2*k*sigma - 44*Q*k*sigma^2 + 11*k*sigma^3 + 36*Q^3 - 54*Q^2*sigma + 32*Q*sigma^2 - 7*sigma^3)/(20*sigma^4);
## where sigma is given in terms of x earlier.
data1 := [seq([lambda1(x), x], x = 0 .. 0.6, 0.1)];

Since lambda1 is NOT defined as a function, but used as one anyway you may wonder why that works (it does!).
To see why, try this version instead:
 

data1 := [seq([lambda1, x], x = 0 .. 0.6, 0.1)];

Notice that you get exactly the same as before because the x inside lambda1 is given the different values of x = 0..0.6 with spacing 0.1.  The first version with lambda1(x) works because any number (e.g. -0.09786536940) will also work as the constant function with that value, thus -0.09786536940(x) = -0.09786536940 for any x.

Since I'm not going to use seq in my code I will turn lambda1 into an actual function of x. That is done by unapply.
After that I can do e.g. lambda1( 0.3) and expect to get what I intended.

2.  N is chosen so that the spacing between the x-values is 10^(-6). The vector V contains all the x-values and is given the datatype float (or explicitly float[8 ] ). Otherwise it would have been datatype=anything. To see that, try:
 

V:=Vector(7,i->0.6/6*(i-1));
op(3,V);

Why do I care about the datatype? Because for speed when V is big (you will have N = 600000) I would like to use hardware float computation (evalhf) and to get the most out of that I avoid data conversion from 'anything' to 'float[8]' by setting datatype=float[8] to begin with.
3. The use of map instead of the easier elementwise operation ~ .
Try

V:=Vector(7,i->0.6/6*(i-1),datatype=float[8]);
sin~(V); # OK
map(sin,V); # OK
evalhf(sin~(V)); # error
evalhf(map(sin,V)); # OK

4. Finally < W | V > creates a matrix with the columns W and V.

Is the vector you give named gln so that gln[i], i = 1..2 are its components?

@tomleslie dsolve, events works in Maple 12, thus also in Maple 13.

Your suggestion [[ gln[1]^2+gln[2]^2, halt]]  should work in principle, but I would like to see the actual system.
Will gln[1]^2+gln[2]^2 ever be found equal to zero, problem being that it is >=0 always?

The solution to that is the obvious: Replace gln[1]^2+gln[2]^2 by gln[1]^2+gln[2]^2-epsilon, where epsilon could be e.g. 1e-7.

There is a dog chasing jogger example in the help page for stop_cond in Maple 8. stop_cond has been superseded by events.
Dog_chasing_jogger_events.mw

For your worksheet to run without error you need kernelopts(floatPi=true).
This option exists in Maple 2015, but is not documented. It doesn't exist in earlier versions. I believe the default is false in 2015.
The option is documented in Maple 2016 and later versions and the default is true.
I found out since I have kernelopts(floatPi=false) in my maple.ini file.