@Suley83 When you copy code like mine and paste it into a worksheet be sure to separate regions where necessary or convenient. Use F3.
I strongly advise you to use 1D math input and worksheet mode. I myself hate 2D-math input.

To change go to Tools/Options/Display/Input Display choose Maple Notation.
Then in the same window go to Interface/Default format for new worksheets. Choose Worksheet.

Finally Press Apply Globally.
New worksheets will now be using Maple Notation (1D).
Here is a minimally revised version of your worksheet:
MaplePrimes20-07-30_SCI_help_Progresses.mw

@Suley83 It is no problem to mimic in Maple what is done in the pdf-file.
So I shall restrict myself to that:
 

RHS; # The right hand sides of odeSys
s1,c1:=selectremove(has,RHS[2],S); # For convenience
### I use the notation in the pdf-file, but will not abuse notation and use F and V in two ways.
### Thus the use of JF and JV:
F:= < s1,0 >;
V:= < eval(c1,z=0), eval(RHS[3],z=0) >;
JF1:=VectorCalculus:-Jacobian(F, [C,Iota]);
JV:=  VectorCalculus:-Jacobian(V, [C,Iota]);
###Assuming that a "disease free equilibrium"  has S = N, C = 0, I = 0  we get for RHS: 
eval(RHS,{S=N,C=0,Iota=0});
###These will all be zero if P1 = P2 = 0 and z = N*v. 
###Under those assumptions we get:
JF:=eval(JF1,S=N);
JF.(-JV)^(-1);
###The eigenvalues for this upper triangular matrix are the values in the diagonal, but will 
### (of course) also be found by: 
LinearAlgebra:-Eigenvalues(%);
R0:=max(%) assuming positive;

Assuming then that z, v, and N are connected through z = N*v we get that
Total(t) = S(t) + C(t) + I(t) is constant if N is interpreted as S(0) + C(0) + I(0).

@Suley83 You could try this:
 

RHS:=evalindets(rhs~([odeSys]),function,f->op(0,f)); # the right hand sides
E:=solve(RHS,{S,C,Iota}); # Equilibria
indets(E,RootOf);
pol:=op([1,1],%); # a polynomial of degree 2 in the unknown _Z
parnames:=indets(RHS,name) minus {C,S,Iota};
parnamesList:=convert(parnames,list);
ics:=S(0)=S0,C(0)=C0,Iota(0)=I0;
params:=[S0,C0,I0,parnamesList[]]; # Must be given numerical values for an illustration
### Using the parameters option to dsolve/numeric:
res:=dsolve({odeSys,ics},numeric, parameters=params);
### Example only. The order is irrelevant if the parameters are given as equations:
Example1:=[S0=100,C0=0,I0=0,k=0.1,o=1,v=.01,w=.5,x=1,z=2,P[1]=.2,P[2]=.3,Q[1]=1,Q[2]=2,N=100];
res(parameters=Example1); # Setting parameters
plots:-odeplot(res,[[t,S(t)],[t,C(t)],[t,Iota(t)]],t=0..500,size=[1200,default],legend=[S,C,Iota]);
eval(E,Example1);
allvalues(%); # Compare to the result below:
res(500);

What is the basic reproduction number you mention?

PS. Notice that I haven't made much of an attempt at "realistic" values of the parameters.

@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.