You didn't succeed in attaching a file. Use the fat green arrow in the editor.

Use PDEtools:-dchange.

For further help upload a worksheet with your actual system.

You say that you tried to attach the worksheet. It didn't work.

Use the big green arrow in the editor for that purpose.

I opened your worksheet in Maple 12. Ran it. Removed the output. Copied all the input and pasted it into Maple 2020.
Ran it then without any problem:
Basically I never save worksheets with output.

@Suley83 I must disappoint you there. I know nothing about stochastic differential equations.

@Suley83 I made a few comments in your worksheet. They are related to the effect of copying from 2D to 1D.
The comments are made in text lines. You get a text line by going to a new input line (insert one above by using Ctrl+k and below by Ctrl+j) and the do Ctrl+t.
So you can make comments in that way, but also by using # in an input line.

SCI_1D_math_with_comments.mw

@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