Your system can be rewritten as a first order system of 4 dependent variables.
I have changed your code some and made extensive use of elementwise operations (using the tilde ~):
 

restart;
A := Matrix(2, 2, [0, -1, 1, -1]);
X := t-> <x(t), y(t)>;
g := X(t)/(t+1)^2;
alpha := -1/4;
X0 := <1, 1>;
B := exp(-2*alpha*(u-t))*Matrix(2, 2, [0, 1, 1, 0]);
eq := diff~(X(t), t) = A.X(t)+int~(B.X(u), u = 0 .. t)+g;
eqs:=expand~(Equate(lhs(eq),rhs(eq)));
S:=[y1(t)=int(exp(u/2)*y(u), u = 0 .. t),x1(t)=int(exp(u/2)*x(u), u = 0 .. t)];
sys1:=subs(rhs~(S)=~lhs~(S),eqs);
sys2:=diff~(S,t);
SYS:={sys1[],sys2[]};
## x(0)=1 and y(0)=1 are given; The last two follow from S.
ics:={x(0)=1,y(0)=1,x1(0)=0,y1(0)=0}; 
res:=dsolve(SYS union ics,numeric);
plots:-odeplot(res,[x(t),y(t)],-0.9..6);

The answer to your actual question about stability I shall leave for you to ponder about.
The plot produced in the code is:

If you look at the answer you get with c>0 you see that e.g. BesselK(0, s/c) appears. For this to make sense c must be different from zero. Thus besides c::real you would need c<>0.  But notice also that BesselK(0,x) is imaginary for x < 0. Thus the answer given for c>0 will probably not be valid for c<0.

That you get the same answer when using c>=0 as when using c>0 is puzzling, though.
But equalities and their negations ( <>) are often ignored by 'assume'.

Try solving for the derivatives and the problem will be obvious:
 

restart;
bvp:={ds*(diff(i(x), x)) = exp(eta(x)), s(x)^3*i(x)*b*r+(1-s(x))^3*a*l*(diff(s(x), x))/s(x)^1.5 = (1-s(x))^3, diff(eta(x), x) = dk*(i(x)-1)/s(x)^p, i(0) = 0, i(1) = 1, s(0) = 1};
sys:=select(has,bvp,diff);
factor(solve(sys,diff~({i,s,eta}(x),x)));

You get for diff(s(x),x) the following:
diff(s(x), x) = s(x)^(3/2)*(s(x)^3*i(x)*b*r+s(x)^3-3.*s(x)^2+3.*s(x)-1.)/(a*l*(s(x)-1.)^3)
As you see the denominator has the factor: (s(x)-1.)^3 which clearly is bad news if you want s(0) = 1.
The numerator is zero at x = 0, since i(0)=0, but it can be written as

s(x)^(3/2)*(s(x)^3*i(x)*b*r+(s(x)-1.)^3)/(a*l*(s(x)-1.)^3);

and again as
s(x)^(3/2)*(s(x)^3*i(x)*b*r/(s(x)-1.)^3 + 1)/al;

Clearly the expression s(x)^3*i(x)*b*r/(s(x)-1.)^3 is very problematic since for this to have no singularity at x = 0 i(x) would have to behave like i(x) = O(x^3) at x = 0.

I don't think that you are doing anything wrong. But there is no solution, which is also clear by inspection of the LinearSolve result.
We see that s[2], s[5], and s[6] must be positive. But then -3*s[2]-s[5]-2*s[6] is negative.
 

solve({x[1]>=0, x[2]>=0, x[3]>=0, x[4]>=0, x[5]>=0, x[6]>=0},[s[1],s[2],s[3], s[4], s[5],s[6]]);

returns [ ].
This could be a bug in LinearMultivariateSystem that it is unable to handle "no solution" in a graceful way.
Try this:
 

restart;
M:=Matrix([[3, 0, 2, 2, 1, 1], [-1, -1, 0, -1, 0, -1], [0, 3, 0, 1, 1, 2], [3, 0, 1, 2, 0, 1], [-1, -1, -1, -1, -1, -1]]);
vec:= Vector[column](5, [3, 1, 0, 0, -2]);
x:=LinearAlgebra:-LinearSolve(M,vec, free = s);
## Debugging SimpleAnd
debug(SolveTools:-Utilities:-SimpleAnd);
SolveTools:-Inequality:-LinearMultivariateSystem({x[1]>=0, x[2]>=0, x[3]>=0, x[4]>=0, x[5]>=0, x[6]>=0},[s[1],s[2],s[3], s[4], s[5],s[6]]);

Notice the green line saying

"enter Utilities:-SimpleAnd, args = -1 = 0, 0 <= -1+2*s[2]+s[5] "

The problem seems to be the false statement -1 = 0 having no parameter.

I have submitted a bug report (SCR).

I will answer your last two questions in my preferred way:
 

restart;
## The right hand sides of your system:
f:=(x,y)->-x+x^2-3*x*y;
g:=(x,y)->-y-5*x*y+y^2;
## The proposed Lyapunov function:
V := (x,y)->x^2/2 + y^2/2;  
## Clearly V(0,0) = 0 and V is positive for (in fact) all (x,y)<>(0,0).
## Now consider:
dV:=<diff(V(x,y),x),diff(V(x,y),y)>.<f(x,y),g(x,y)> assuming real;
## For local stability at (0,0) we need only consider the lowest order terms of dV.
## Those are of order 2:
mtaylor(dV,[x,y],3); # -x^2-y^2: Clearly negative for (x,y)<>(0,0)

Conclusion: dV is negative in a neighborhood U of (0,0) if (x,y) <>(0,0).
Thus (0,0) is (locally) asymptotically stable.
To give an explicit neighborhood we can make a few estimates:
 

dV2:=mtaylor(dV,[x,y],3);
dV3:=expand(dV-dV2);
dV3a:=map(abs,dV3);
## Crude, but correct:
dV3a<=(abs(x)+3*abs(y)+5*abs(x)+abs(y))*(x^2+y^2);

Thus if
6*abs(x)+4*abs(y) <= k and k < 1;
 we have dV = dV3+dV2 <= dV3a+dV2 <=
(6*abs(x)+4*abs(y))*(x^2+y^2)+dV2 = -(1-k)*(x^2+y^2) , so dV < 0 if also (x,y)<>(0,0).

So a neighborhood U that works is defined by the inequality: 6*abs(x)+4*abs(y) <= k, where k < 1.

The basin of attraction of (0,0) is bounded by the cyan colored separatrices to the 3 saddle points in this phase plane image:

Since you have

X:=x(tau): 

you cannot use X in your IC. You must use x.
Furthermore w1 needs to be given a value for the numerical solution to work:
 

IC:=x(0)=0,D(x)(0)=0;
##
sys1:={sys,IC}:
dsolve(eval(sys1,w1=1),numeric,method=rkf45); # Example

Whether a symbolic solution ever appears you will see if you are more patient than I am.

You can check your results with dsolve/numeric.
If you remove the method option method=classical[rk4], then you get the default method=rkf45 (then also remove the stepsize option). Using the parameters option allows you to change the values of x0, k, and f in a very simple manner as below.

## Assuming that m, F, and Omega have been defined as in your worksheet.
k:='k': f:='f':
### Using method=classical[rk4] for direct comparison
solRK4:=dsolve({ode,x(0)=x0,D(x)0)=0},numeric,method=classical[rk4],stepsize=0.001,parameters=[x0,k,f]);
solRK4(parameters=[3,0.1,0]); #Setting parameters
plots:-odeplot(solRK4,[x(t),diff(x(t),t)],0..10);
plots:-odeplot(solRK4,[t,1/2*diff(x(t),t)^2+V(x(t))],0..10); #Energy plot
plots:-odeplot(solRK4,[[t,x(t)],[t,diff(x(t),t)]],0..10,size=[800,default]);

You can set parameters as done above as you wish.

You may try depends.

depends(expr,a);

where expr is the polynomial or equation.

There could very likely be a singularity, not only a numerical problem.
I made the following changes:
Digits:=15:
# In the loop used this try..catch structure:
 

try
   sh0:=SH(0):
   sv0:=SV(0):
   Lx0:=LA(0):
   print('i'=i);
 catch:
   print("Breaking",'i'=i); break
 end try;

around the 3 assignments.
Furthermore, in the dsolve command I put in stiff=true, maxfun=10^6, and removed the range argument.
Then the loop ran just once. The break came at i = 2 probably because of an actual singularity strongly indicated here:

sol(parameters);
odeplot(sol,[x,L(x)],4.238507..H-He);
odeplot(sol,[x,sh(x)],4.238507..H-He);
odeplot(sol,[x,sv(x)],4.238507..H-He);
odeplot(sol,[x,u(x)],4.238507..H-He);

Looking at the output from the procedure definition of S1 (the blue stuff) I see psi inside the procedure body although the input appears to have phi (actually varphi).
My guess is that you had a first version defined like that. Then you realized that that was wrong and changed the input. You didn't execute the revised definition, however.

Just remove 'numeric':
 

int(piecewise(0 < -k5*x1*x2+k1*x1 and -k5*x1*x2+k1*x1 < 1, 1, -k5*x1*x2+k1*x1 <= 0 and 1 <= -k5*x1*x2+k1*x1, 0)*piecewise(0 < -k5*x1*x2+k2*x2 and -k5*x1*x2+k2*x2 < 1, 1, -k5*x1*x2+k2*x2 <= 0 and 1 <= -k5*x1*x2+k2*x2, 0), k1 = 0 .. 1, k2 = 0 .. 1, k5 = 0 .. 1, x1 = 0 .. 1, x2 = 0 .. 1);
## Answer 7/12
int(piecewise(0 < x+y and x+y < 1, 1, x+y <= 0 and 1 <= x+y, 0)*piecewise(0 < y*x^2 and y*x^2 < 1, 1, y*x^2 <= 0 and 1 <= y*x^2, 0), x = 0 .. 1, y = 0 .. 1);
## Answer 1/2

That numeric integration takes a long time or has its problems is likely due to your integrands not being continuous.

Maybe it is a slight exaggeration to call your problem nonsense.
But my reason is this. Your system can in vector form be written as
X'(t) = f(t, X(t)), X(t0) = X0  (*)
where X and X0 are vectors, in your case of dimension 2, and f is a vector valued function.
The initial value problem (*) is equivalent to the integral equation
X(t) = X0 + int( f(s, X(s) ), s=t0 .. t);
As the value of an integral stays the same if the integrand is changed on a set of measure zero, your problem is equivalent to the simple problem
 

sys:={diff(x(t),t)=y(t), diff(y(t),t)=-x(t), x(0)=0,y(0)=1};

In the definition of c a final parenthesis ')'  is missing. After that no complaint from the parser.

Isn't it a problem that your function has infinitely many singularities on the positive real axis
(s[n], n=1..infinity) where s[n] -> infinity as n - > infinity
?
To see that there are singularities as described (and as expected because of the tan function) you could do:
 

restart;
F:=1/(s^2*(1+(1-1/(2+2*s))/s)*(1+tan((1/2)*Pi/sqrt((1-1/(2+2*s))/s))/sqrt((1-1/(2+2*s))/s)));
F2:=normal(convert(F,sincos));
A:=op([1,1],indets(F2,specfunc(sin)))/Pi;
limit(A,s=infinity);
sp1,sp2:=solve(A=p,s);
evalf(seq(sp1,p=1..10));
dF2:=denom(F2);
seq(fsolve(dF2=0, s=eval(sp1,p=n)..eval(sp1,p=n+1)),n=0..9);

The sine term in the denominator dominates the cosine term when s > 1 and tends to zero as s->infinity. Thus the singularities will be close to those values of s that make the argument to sine an integral multiple of Pi. Those values are given by sp1 for p = 1..infinity.
The first 10 actual singularities are found in the last command by fsolve.

Thanks to acer for suggesting major parts of the following approach.

Making use of dsolve 20000 times takes a while. So if the system could be turned into an initial value problem we can use the parameters option.

Although your system is really a boundary value problem since it has four functions with conditions given at H-He and one function (u) with u(0) = 0 we can consider instead posing the condition u(H-He)=0.
Your system is special. The ode for u has the form diff(u(x),x) = f( Jir(x) ) = g(x), i.e. it doesn't depend on u(x).
Thus the u(x) is just obtained by integration, u(x) = int(g(xi), xi = H-He .. x) (not that we shall do that literally).
The function u we actually want (ua, say) is then given by ua(x) = u(x)-u(0).
I have modified other parts, but that is really not important.
 

restart;
Digits:=15:
phi0:=0.72:
rho0:=646.8:
E:=100*10^6:
nu:=1/3:
pc0:=0.3*10^6:
pv0:=0.6*10^6:
g:=9.80665:
Mp:=rho0*100/(10^6*365*24*60*60):
lambda:=E*nu/((1+nu)*(1-2*nu)):
mu:=E/(2*(1+nu)):
K:=lambda+2/3*mu:
eta:=100*10^6*10^6*365*24*60*60:
a:=1.15:
A:=-2/3+a*sqrt(3)/3:
B:=-1/3-a*sqrt(3)/3:
m:=5:
n:=3:
Ts:=60*10^6*365*24*60*60:
Tf:=2*Ts:
####
Md:=Mp*t:
phi:=1-(1-phi0)/Jir(x);
pc:=pc0*(ln(phi)/ln(phi0))^m;
mp:=-(m*(1-phi0)*pc)/(Jir(x)*phi*ln(phi));
Gp:=-(K+2*sqrt(3)/3*mu*a)/(K+a^2*mu+mp);
beta11:=K-2/3*mu;
beta21:=K+4/3*mu;
beta12:=(K-sqrt(3)/3*mu*a)*Gp+beta11;
beta22:=(K+2*sqrt(3)/3*mu*a)*Gp+beta21;
####
Delta1:=1/(rho0*g);
H:=beta21*Delta1*(1-exp(-g/beta21*Md));
Lambda:=1-rho0*g/beta21*(H-x);
####
sh:=beta11*ln(Lambda);
sv:=beta21*ln(Lambda);
Le:=evalf(exp(-pc0/(K+2*sqrt(3)/3*mu*a)));
Jire:=1;
He:=beta21*Delta1*(1-Le);
Te:=fsolve(H=He,t);
she:=eval(sh,{x=0,t=Te});
sve:=eval(sv,{x=0,t=Te});
ue:=-Mp*g/beta21*He;
unassign('Lambda','sh','sv','H'):
###############################
sys_ode:=diff(sv(x),x)=rho0*g/Lambda(x),diff(sh(x),x)=beta12/beta22*rho0*g/Lambda(x),diff(Lambda(x),x)=rho0*g/beta22,diff(u(x),x)=-Mp*g/beta22,diff(Jir(x),x)/Jir(x)=-Gp*diff(Lambda(x),x)/Lambda(x);
### Your ics has u(0)=0. The new ones:
icsNew:={sv(H-He)=sve,sh(H-He)=she,Lambda(H-He)=Le,Jir(H-He)=Jire,u(H-He)=0};
resNew:=dsolve({sys_ode} union icsNew,numeric,[sv(x),sh(x),Lambda(x),Jir(x),u(x)],parameters=[H],
        output=listprocedure,abserr=1e-14,relerr=1e-12); #optimize may help a little
SV,SH,JIR,U:=op(subs(resNew,[sv(x),sh(x),Jir(x),u(x)]));
resNew(parameters=[He]); #Setting the parameter H to He for a test:
SV(0),SH(0),JIR(0),U(0);
#######
nn:=20000:
CC1:=Matrix(datatype=float[8]):
CC2:=Matrix(datatype=float[8]):
TR1:=Matrix(datatype=float[8]):
TR2:=Matrix(datatype=float[8]):
EN1:=Matrix(datatype=float[8]):
EN2:=Matrix(datatype=float[8]):
####
t0:=time():
dt:=(Ts-Te)/nn:
t:=Te:
HH:=He:
for i from 1 to nn+1 do
 resNew(parameters=[HH]); #Setting the parameter H to HH
 CC2(i,1..2):=Array([t/(10^6*365*24*60*60),HH]):
 sh0:=SH(0):
 sv0:=SV(0):
 TR2(i,1..2):=Array([-1/3*sh0-2/3*sv0,evalf(sqrt(3)/3*(sh0-sv0))]):
 Jir0:=JIR(0):
 phix0:=1-(1-phi0)/Jir0:
 EN2(i,1):=pc0*(ln(phix0)/ln(phi0))^m:
 up:=U(HH-He)-U(0); # Since U is just an integral.
 uh:=ue+up:
 fvp:=evalf[30](A*sh0+B*sv0-(pv0*(ln(phix0)/ln(phi0))^n)):
 if fvp>=0 then
  break
 else
  dH:=dt*(Mp/rho0+uh):
  HH:=HH+dH:
  t:=t+dt:
 end if:
end do:

time()-t0;
i;

The time used on my computer was 3.140s and i evaluates to 742.
MaplePrimes18-08-21_odesys_inits.mw