@mmcdara I believe that your ode is unstable in the forward direction, i.e. as t -> infinity, but stable in the opposite direction.
If I'm right and you are interested in a solution on the finite interval t = 0..t0, then start the solution at t0.
Since you only know X and D(X) at 0 you could try use a shooting method.
But even if you succeed in find starting values (faking not knowing U) for a solution starting at t = t0 (t0 > 0), then that solution will only be good on t = 0..t0. Just try using X(t0) = U(t0) and D(X)(t0) = D(U)(t0) as your initial contions.
You get perfect results for t = 0..t0, but for t > t0 not so.
An example:
Here I chose t0 = 1 and plottet X(t) on 0..1.35:

@pik1432 It is easy to write down several ode systems all with equilibria coming from solving your algebraic equations.
Just multiply say 3 of these equations by minus 1 and the equilbria are the same, but the jacobian and its eigenvalues are not.
Thus stability here must refer to the solving method and obviously must refer to a numerical method as used by fsolve, not an exact as used by solve.
##################################################################################
A trivial one dimensional example is
diff(x(t),t) = x(t) which has jacobian 1 at the equilibrium x = 0 and eigenvalue 1: 0 is unstable.
diff(x(t),t) = -x(t) which has jacobian -1 at the equilibrium x = 0 and eigenvalue -1: 0 is stable.
Yet in both cases determining the equilibrium consists in solving the same basic equation x=0.
 

@pik1432 I forgot that I used expand on your f9n10 as I see that also acer does.

Apparently you are using some hidden code. That doesn't exactly help.

Are you using 2D input?
If so then don't.
Go to Tools/Options/Display/Input display.  Pick Maple Notation (i.e 1D).
Then still with the Options window open go to Interface/Default format for new worksheets. Pick Worksheet.
Now press the button Apply Globally.

Open a new worksheet. I hope you like it better. I do!

If you use solve instead you will get a solution.

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.