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