restart:

with(plots):

eq1 := diff(f(x), x, x, x)+(1/2)*cos(alpha)*x*(diff(f(x), x, x))+(1/2)*sin(alpha)*f(x)*(diff(f(x), x, x)) = 0;

eq2 := diff(g(x), x, x)+diff(g(x), x)+(diff(g(x), x))*(diff(h(x), x))+cos(alpha)*x*(diff(g(x), x))+sin(alpha)*f(x)*g(x) = 0;

eq3 := diff(g(x), x, x)+diff(h(x), x, x)+1/2*(cos(alpha)*x+sin(alpha)*f(x)) = 0;

ics := f(0)=0, D(f)(0)=1, (D@@2)(f)(0)=a[1], g(0)=1, D(g)(0)=a[2], h(0)=1, D(h)(0)=a[3];

#bcs:=f(x) , g(x), h(x) tends to 0 as x tends to infinity

sys := {eq1, eq2, eq3, ics}:

# dsolve(sys); # It's unlikely that a formal solution does exist

Warning,  computation interrupted

params := convert(indets(sys, name) minus {x}, list)

sol := dsolve(sys, numeric, parameters=params);

# example of use:
#
# Step 1: instanciate parameters

data := [$1..4]:

sol(parameters=data);

# Step 2 (illustration): plot some quantity

display(
  odeplot(sol, [x, f(x)], x=0..1, color=red, legend=typeset('f(x)'))
  , odeplot(sol, [x, g(x)], x=0..1, color=green, legend=typeset('g(x)'))
  , odeplot(sol, [x, h(x)], x=0..1, color=blue, legend=typeset('h(x)'))
)

# Method 1
#
# Assuming you want to replace the ICs which contain the higher derivative order,
# and assuming that L is a reasonnable approximation of +oo

L := 10:
`ICs+BCs` :=  f(0)=0, D(f)(0)=1, g(0)=1, h(0)=1, f(L)=0, g(L)=0, h(L)=0;

# REMARK:
#
# As soon as your differential susyemen containsBCs, the use of option "parameters"
# indsolve/numerics is forbidden.
# Thus you have to instanciate the parameters BEFORE using dsolve.
SYS := {eq1, eq2, eq3, `ICs+BCs`}:

SOL := dsolve(eval(SYS, params=~data), numeric)

# Plots

display(
  odeplot(SOL, [x, f(x)], x=0..L, color=red, legend=typeset('f(x)'))
  , odeplot(SOL, [x, g(x)], x=0..L, color=green, legend=typeset('g(x)'))
  , odeplot(SOL, [x, h(x)], x=0..L, color=blue, legend=typeset('h(x)'))
)

# REMARK:
#
# This naive method doesn't always work


L := 100:
`ICs+BCs` :=  f(0)=0, D(f)(0)=1, g(0)=1, h(0)=1, f(L)=0, g(L)=0, h(L)=0;

SYS := {eq1, eq2, eq3, `ICs+BCs`}:
SOL := dsolve(eval(SYS, params=~data), numeric, maxmesh=10^4)

Error, (in dsolve/numeric/bvp) Newton iteration is not converging

# Method 2
#
# Re-parameterize the problem by using a 1-to-1 map from (0, +oo) into (0, 1).
#
# Here I use the change phi := x -> 1/(1-x) in such a way that +oo is mapped to 1/

edo := eval(eqs, {f(x)=f(phi(x)), g(x)=g(phi(x)), h(x)=h(phi(x))}):

convert(select(has, indets(edo, function), D), list):
Changes := % =~ subsop~(-1=z, %):
eval(eval(edo, Changes), {f(phi(x))=f(z), g(phi(x))=g(z), h(phi(x))=h(z)}):

EDO := (numer@lhs@simplify)~(convert(eval(eval(%, phi=(x->1/(1-x))), x=z), diff));

`ICs+BCs` :=  f(0)=0, D(f)(0)=1, g(0)=1, h(0)=1, f(1)=0, g(1)=0, h(1)=0;

SYS := {EDO[], `ICs+BCs`}:
SOL := dsolve(eval(SYS, params=~data), numeric, maxmesh=10^4)

display(
  odeplot(SOL, [1/(1-z), f(z)], z=0..0.999, axis[1]=[mode=log], color=red, legend=typeset('f(x)'), labels=[x, ""])
  , odeplot(SOL, [1/(1-z), g(z)], z=0..0.999, axis[1]=[mode=log], color=green, legend=typeset('g(x)'))
  , odeplot(SOL, [1/(1-z), h(z)], z=0..0.999, axis[1]=[mode=log], color=blue, legend=typeset('h(x)'))
)

plot(
  (subs(M = 1.3015, A)-subs(M = 1.2431, A))/subs(M = 1.2431, A)*100
  , x = -5 .. 0
  , axes = boxed
  , colour = red
  , tickmarks=[default, [seq(k=sprintf("%d%%", k), k in [seq](-60..0, 10))]]
  , title="(A-Aref)/Aref*100%"
  , gridlines=true
)

restart:

with(Statistics):

B := RandomVariable('Beta' (alpha, beta))

BSD := x__min + (x__max-x__min)*B

# PDF of BSD

f := PDF(BSD, x)

# PDF of BSD in the form given in your reference

eval(f) assuming x > x__min, x__max > x__min, x < x__max:

subs(x__max-x__min=h, %):
map(simplify, %) assuming alpha > 0, beta > 0, h > 0:
simplify(%):

f__BSD := subs(h=x__max-x__min, %);
 

# A closed form of the CDF does exist contraryto what is said in your reference

F := CDF(BSD, x);
print():

# After accounting for assumptions:

eval(F)  assuming x > x__min, x__max > x__min, x < x__max
 

# Direct way (I don't know how to simplify this expression)


'mean(BSD)' = simplify(simplify(Mean(BSD), GAMMA), GAMMA);

# alternative way (using the linearity of the "Expectation" operator

'mean(BSD)' = x__min + (x__max-x__min)*Mean(B);

# Direct way (I don't know how to simplify this expression)

'variance(BSD)' =  simplify(Variance(BSD), GAMMA):

# alternative way (using the linearity of the "Expectation" operator

'variance(BSD)' = (x__max-x__min)^2*Variance(B);

# Direct way (fails)

'mode(BSD)' = simplify(Mode(BSD), GAMMA);

# alternative way (obvious result)

'mode(BSD)' = x__min + (x__max-x__min)*Mode(B);

Warning, solutions may have been lost

B   := RandomVariable('Beta' (alpha, beta)):
BSD := x__min + (x__max-x__min)*B:

# f and F above are used to be copied-pasted into SubjectiveBeta

f   := PDF(BSD, x);
F   := CDF(BSD, x) assuming x__max > x__min

SubjectiveBeta := proc(A, B, P, Q)
   Distribution(
     PDF=unapply(
       eval(
         piecewise(
           (x-x__min)/(x__max-x__min) < 0, 0
           , (x-x__min)/(x__max-x__min) < 1
           , ((x-x__min)/(x__max-x__min))^(-1+alpha)*(1-(x-x__min)/(x__max-x__min))^(-1+beta)/Beta(alpha, beta)
           , 0
         )/(x__max-x__min)
         , [alpha=A, beta=B, x__min=P, x__max=Q]
       )
       , t
     ),
     CDF=unapply(
       eval(
         piecewise(
           (x-x__min)/(x__max-x__min) < 0, 0
           , (x-x__min)/(x__max-x__min) < 1,
               ((x-x__min)/(x__max-x__min))^alpha
               *
               hypergeom([alpha, 1-beta], [alpha+1], (x-x__min)/(x__max-x__min))/(Beta(alpha, beta)*alpha)
           , 1
         )
         , [alpha=A, beta=B, x__min=P, x__max=Q]
       )
       , t
      ),

      Mean = eval(
               x__min+(x__max-x__min)*alpha/(alpha+beta)
               , [alpha=A, beta=B, x__min=P, x__max=Q]
             ),

      Variance = eval(
                   (x__max-x__min)^2*alpha*beta/((alpha+beta)^2*(alpha+beta+1))
                   , [alpha=A, beta=B, x__min=P, x__max=Q]
                 ),

      Mode = eval(
               x__min+(x__max-x__min) * Mode('Beta'(A, B))
               , [alpha=A, beta=B, x__min=P, x__max=Q]
             ),

      Conditions = [A > 0, B > 0, Q > P],

      RandomSample = proc(N::nonnegint)
                       P +~ (Q-P) *~ Sample('Beta'(A, B), N)
                     end proc
   )
end proc:

BSD := (alpha, beta, xmin, xmax) -> RandomVariable(SubjectiveBeta(alpha, beta, xmin, xmax)):

SubjBeta := BSD(alpha, beta, x__min, x__max);

PDF(SubjBeta, x);

CDF(SubjBeta, x);

Mean(SubjBeta);
Variance(SubjBeta);
Mode(SubjBeta);

SubjBetaNum := BSD(3, 1, -1, 4);

PDF(SubjBetaNum, x);

CDF(SubjBetaNum, x);

Mean(SubjBetaNum);

Variance(SubjBetaNum);

Mode(SubjBetaNum);

Histogram( Sample(SubjBetaNum, 1000) )

kernelopts(version);

# What you want to "dsolve" is:

ToDsolve := [ subs(a1,eq1), bcs, subs(a1,eq2), bcs1 ]:

print~(ToDsolve):
 

# ToDsolve contains some formal quantities:

FQ := convert(indets(ToDsolve, name) minus {f, theta, x}, list)

# Give them numerical values and look to what happens

dsolve(eval(ToDsolve, FQ=~1), numeric)

Error, (in dsolve/numeric/bvp/convertsys) too many boundary conditions: expected 5, got 6

# The message means that you have 3 bcs for theta as theta verifies a 2nd order ode
# It's upt to you to fix that.
# I just show you that once fixed thete is no longer any error:
#
# Attempt 1

attempt := subsop(7=NULL, ToDsolve):
print~(attempt ):
dsolve(eval(attempt , FQ=~1), numeric)

# Attempt 2

attempt := subsop(8=NULL, ToDsolve):
print~(attempt ):
dsolve(eval(attempt , FQ=~1), numeric)