If you do

f1:=1+x;
g:=x->f1;
g(8);
## then you see that g(8) returns 1+x. In fact any input to g would return 1+x.

Executing a procedure definition as in g:=x->f1; does not cause evaluation of the body (in this case just f1).
Evaluation doesn't take place before the procedure is called (applied), as in g(8).
The first thing that happens when you do g(8) is that the x's appearing in the body (here f1) are replaced by 8.
Since there hasn't been any evaluation yet, there won't be any x to be replaced.
After replacement, evaluation occurs. f1 evaluates to 1+x and that is what is returned.

That is why you need unapply in your case.

This example may illustrate this further:
h:=x->f1*x;
h(8);
#There was an x to be replaced by 8, so it was. Then evaluation and you get 8+8*x.
Notice also that just g(); results in 1+x, but h(); results in an error, because there is an actual x to be replaced inside h.

Since you have a symbolic S you need to give the dependent variable:

dsolve({eqn0,bcs0},f0(y));

# Answer:    f0(y) = -(1/2)*y^3/S^3+(3/2)*y/S

This would be unnecessary had S had a numeric value.

The syntax for symbolic solution of a system of pdes with ibcs is different from the syntax for numerical solution.
For a symbolic solution you give the pdes and the ibcs together as the first argument to pdsolve (as a set) as you would do for odes.
For a numerical solution you give the system as the first argument and the ibcs as the second argument and you need to add the argument 'numeric'.
There is a numeric example of a system in the help page for pdsolve,numeric at the very end.
It seems that you are right there is no example of solving a system symbolically with conditions.
## Trying to construct a simple example:

pde:=diff(u(x,t),t,t)=diff(u(x,t),x,x); #The wave equation
pdsolve(pde);
sys:=diff(u(x,t),t)=w(x,t),diff(w(x,t),t)=diff(u(x,t),x,x); #Corresponding system
pdsolve({sys}); #OK so far
pdsolve({pde,u(x,0)=sin(x),D[2](u)(x,0)=1}); #Looks good
diff(%,t); # so this would be w(x,t), however:
pdsolve({sys,u(x,0)=sin(x),w(x,0)=1}); #Cannot do it!

You should not use 'with' inside a procedure.
In the help page for 'with' ( ?with ) you find an explanation why and the example

p := proc()
    with( ListTools );
    Reverse( [ 1, 2, 3 ] )
end proc:
p();

A solution is to use 'uses':

p := proc()
    uses ListTools;
    Reverse( [ 1, 2, 3 ] )
end proc;  #semicolon so you see the effect 'uses' has when the procedure definition is executed
p();

You need to use unapply.
Furthermore it is better to use eval than subs, but not essential. Your "weird" substitution 0=0 only has the effect of causing an evaluation, which is an essential difference between subs an eval. subs doesn't evaluate after substitution.

restart;
changering := proc(Equation1, f3,g3) local g1,f1,h; global x,y,f,g;
g1 := unapply(f3,x,y);
f1 := unapply(g3,x,y);
eval(Equation1,{g=g1,f=f1})
end proc:
Eq1 := f(x,g(x,y)) + f(x,y);
h2 := changering(Eq1,x+y, x+y);
g1 := (x,y)-> x+y;
f1 := (x,y)-> x+y;
h:=subs(g=g1, Eq1);
h:=subs(f=f1, h);
%;
##Compare to:
eval(Eq1,{g=g1,f=f1});

We can certainly conclude that the entirely symbolic result 'sc', i.e. the result obtained from solve before values of parameters are inserted, does not represent correctly all 3 roots for all values of theta, R, h, H.
In fact it is entirely incorrect for your choice of the parameters 'va': All 3 roots are the same and wrong.
Even when inserting your values for R, H, and h, but not theta the solve result is totally wrong for theta=Pi/6.

It is in general safer to insert values (like va) before solving using solve. The (again) symbolic result obtained can then be compared to solving numerically using fsolve with the option 'complex'.

restart;
vol := -R^2*Pi*H*(-3*R*C^2*sin(theta)-3*H^2*cos(theta)^2*sin(theta)*R+6*R*C*sin(theta)*cos(theta)*H-sin(theta)*R^3+sin(theta)*R^3*cos(theta)^2+C^3-3*C^2*cos(theta)*H-cos(theta)^3*H^3+3*H^2*C*cos(theta)^2+3*R^2*H*cos(theta)^3+3*R^2*C-3*R^2*H*cos(theta)-3*cos(theta)^2*R^2*C)/(3*(H^2*cos(theta)^2-R^2+cos(theta)^2*R^2)^(3/2));
V := (1/3)*Pi*R^2*H;
vp := Pi*h*R^2*(3*H^2-3*H*h+h^2)/(3*H^2);
eq := V-vp = vol;
va := theta = (1/6)*Pi, R = 2, H = 20, h = 14;
##Up til now I was just repeating your code####
eq_va:=eval(eq,{va}); #The concrete equation in the unknown C
solve(eq_va,C); #Symbolic result
evalf(%);
fsolve(eq_va,C,complex); #Solving numerically
## The symbolic result is confirmed by the numerical
##Now the converse: insertion of parameters last
#Parameters not inserted:
sc := solve(eq,C): #sequence now
eval({sc},{va});
simplify(%); #All 3 equal
evalf(%); # ... and wrong
###########Partial insertion
eval({sc},{va[1..3]}); #All but h inserted
simplify(%) assuming h<20; #WRONG result (unfortunately, because h = 14 < 20)
simplify(%%) assuming h>20; #Correct at least for h = 14 !
eval(%,va[4]);
evalf(%);
###Partial insertion: all but theta. Fails entirely in case theta = Pi/6
EQ:=eval(eq,{h=14,H=20,R=2});
res:=solve(EQ,C);
indets({res},RootOf); #This Rootof is just a square root the same at all occurrences
allvalues(%);
eval({res},theta=Pi/6);
simplify(%); #Wrong

If you think that this linear homogeneous system with homogeneous boundary conditions has a nontrivial solution, then you could try finding one by using an approximate solution, which basically just has to be nonzero, as in

res:=dsolve({sys,cond},numeric,approxsoln=[U(xi)=xi,y(xi)=(1-xi)*(xi-0.8)]);

(I used xi instead of ksi: it looks better in output).

I do get consistent shapes, but the sizes vary dependent on the setting of Digits, which is troubling.

Since the problem is linear and homogeneous certainly if (U,y) is a (nontrivial) solution so is any constant multiple, so the fact that the maximal values of the solution is very small doesn't prove that it is wrong. However, it would be nice to get a solution where the maximal value for y is 1 instead of 4*10^(-42).

Your system may be 'close' to a system that actually has nontrivial solutions. You could try introducing a parameter L and then add a nonzero boundary condition for U. You could hope for nontrivial solutions (y,U) for L close to 1. Indeed that turns out to be the case.

restart;
Digits:=15:
sys := diff(U(xi ), xi ) = (y(xi )-(1.5*(1+0.01*xi ))*U(xi )/xi )/(1.5+0.015*xi +0.002*xi ^2),
diff(y(xi ), xi ) = U(xi )*(0.002+(1.5*(1+0.01*xi ))*(0.002*xi ^2)/(xi ^2*(1.5+0.015*xi +0.002*xi ^2))-19.3^2)-y(xi )*(0.002*xi ^2)/(xi *(1.5+0.015*xi +0.002*xi ^2));

cond:= y(0.8) = 0, y(1) = 0;
res:=dsolve({sys,cond},numeric,approxsoln=[U(xi)=xi,y(xi)=(1-xi)*(xi-0.8)]);
res(.95);sys}
plots:-odeplot(res,[[xi,y(xi)],[xi,U(xi)]]);
###Now introduce a parameter L
sysL:=subs(y(xi)=L*y(xi),U(xi)=L*U(xi),convert({sys},D));
resL:=dsolve(sysL union {cond,U(1)=1},numeric,approxsoln=[L=1,U(xi)=xi,y(xi)=(1-xi)*(xi-0.8)]);
resL(.9); #Notice that L = 1.00118260702656 i.e. close to 1.
plots:-odeplot(resL,[[xi,y(xi)],[xi,U(xi)]]);

If I understand this correctly, there is no real problem.
Example (silly):

p:=proc(x::numeric) local q,k,r;
  q:=proc(s) s^2 end proc:
  r:=x;
  for k from 1 to 5 do
    r:=q(r)
  end do;
  evalf(sin(r))
end proc;
p(2);
p(89);

You had a problem with parentheses in the assignment to b. I changed it to a version which may not be what you intended. Also the ´global` declaration should be inside the procedure.
It is OK with the explicit return statement at the end, but being at the end it is superfluous.

restart;

a:=0.081819221;
PI:=3.1415926535897932384626433832795;
Ecce:=proc(lt) local lat,b,t; global a,PI;
  lat:=lt*(PI/180);
  b:=(1-a*sin(lat))/(1+a*sin(lat))^(a/2);
  t:=ln((tan(PI/4)+lat/2));
  3437.7468*t*b
end proc;
Ecce(45.2112);

You will very easily run into a nonuniqueness problem when using fsolve on your system.
To examine it closer notice that you can express yp2 in terms of yp4. You can find an equation for yp4 expressed in terms of Y[i], i=1..4.
If you do that you get a version of p as follows.
The globals defined were used for diagnostical purposes.

gg:=1:
eps:=2.:
epsmax:=10:
p2:=proc(N,t,Y,YP)
  local eq,yp4; global gg,eps,gg1;
  YP[1]:=Y[2];
  YP[3]:=Y[4];
  eq:=Y[2]*yp4*sin(Y[1]*Y[3])-5*(yp4*Y[2]*sin(Y[1]^2)+cos(yp4*Y[3])-sin(t))*Y[4]*cos(Y[1]^2)/Y[3]+t^2*Y[1]*Y[3]^2 = exp(-Y[3]^2);
  gg:= fsolve(eq,yp4=gg-eps..gg+eps);
  if has(gg,fsolve) then gg:=gg1; eps:=min(eps+1,epsmax) end if;
  gg1:=gg;
  YP[4]:=gg;
  YP[2]:= -(YP[4]*Y[2]*sin(Y[1]^2)+cos(YP[4]*Y[3])-sin(t))/Y[3];
  0
end proc:
res2:=dsolve(numeric,procedure=p2,initial=Array([1,1,2,2]),number=4,procvars=[x1(t),diff(x1(t),t),x2(t),diff(x2(t),t)],start=0);
res2(0.28);
gg; # 916.---
eps; # 10
You have the problem that the equation eq, which has the form a*z+b*cos(d*z)=f might easily have several solutions for z. You have to pick the right one, and if it happens that the previously found root is a double root then the this double root may become two roots or none, so what to do?

As a small illustration try the plot
plot(z+20*cos(z)-1,z=-30..40);

or the animation
plots:-animate(plot,[a*z+20*cos(z)-1,z=-30..30],a=.9..1.1,size=[1800,400]);

restart;
A1 := Matrix([[a11, a12, a13], [a12, a22, a23], [a13, a23, a33]]);
 A2 := Matrix([[A], [B], [C]]);
 A3 := Matrix([[15], [0], [0]]);
eq := A1.A2=A3; #Your equation
 
LinearAlgebra:-LinearSolve(A1,A3);

Try this:

restart;
n := 3;
y[0] := x-> 1+x+(1/2)*x^2;
for m from 0 to n do
  y[m+1]:=unapply( y[m](x)+int(-(1/2)*(s-x)^2*(diff(y[m](s),s,s,s)-y[m](s)), s = 0 .. x),x)
end do;
ode:=diff(Y(x),x,x,x)-Y(x)=0;
dsolve({ode,Y(0)=1,D(Y)(0)=1,(D@@2)(Y)(0)=1});
sort(collect(y[n](x),x),x,ascending);
mtaylor(exp(x),x=0,11);
%%-%;

In early versions of Maple all packages were tables. That is true of plots still in Maple 8.

So replace plots:-pointplot with plots[pointplot] . That can even be done in later versions.

Before getting to the parameter you are talking about (actually your worksheet doesn't contain any) you should address the problem that you have: an error message from `dsolve/numeric/bvp/convertsys`.

A solution to this problem is to differentiate the equation not containing the 4th derivatives. In addition use the undifferentiated version to extract an extra boundary condition, because you only have 10, you need 4+4+3 = 11.

restart;
##I should excuse the unnecessarily complicated look of dsys3. I converted to 1-D math input and didn't bother to clean up the mess:
dsys3 := {8*(diff(f2(x), x, x, x, x))+9*(diff(f2(x), x, x))+10*f2(x)+11*(diff(f1(x), x, x, x))+12*(diff(f1(x), x))+13*(diff(f3(x), x, x))+14*f3(x)+f3(x)*f3(x)+(diff(f3(x), x))*(diff(f3(x), x))+(diff(f3(x), x, x))*f3(x) = 0, 16*(diff(f3(x), x, x, x, x))+18*(diff(f3(x), x, x))+19*(diff(f3(x), x, x))+22*(diff(f1(x), x))+23*(diff(f1(x), x))+24*(diff(f2(x), x, x))+25*f2(x)+26*f2(x)+27*f3(x)+29*f3(x) = 0, 2*(diff(f1(x), x, x))+3*(diff(f2(x), x, x, x))+4*(diff(f2(x), x))+6*(diff(f3(x), x))+7*f1(x)+(diff(f3(x), x, x))*(diff(f3(x), x))+(diff(f3(x), x))*f3(x)+(diff(f3(x), x))*f3(x) = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, ((D@@1)(f2))(0) = 0, ((D@@1)(f2))(1) = 0, ((D@@1)(f3))(0) = 0, ((D@@1)(f3))(1) = 0};
dsys3;
sys,bcs:=selectremove(has,dsys3,diff);
solve(sys,{diff(f3(x),x$4),diff(f2(x),x$4),diff(f1(x),x$3)});
sys[1];
sys[2];#In my case the one I want to differentiate
sys[3];
sys2:=diff(sys[2],x);
#Check that we can solve for the highest derivatives (and defining the new system:
SYS:=solve({sys[1],sys2,sys[3]},{diff(f3(x),x$4),diff(f2(x),x$4),diff(f1(x),x$3)});
bcs;
#Finding one more boundary condition:
convert(sys[2],D);
eval(%,x=0);
bcs1:=eval(%,bcs);
## Now we can solve the problem:
dsol:= dsolve(SYS union bcs union {bcs1}, numeric);
plots:-odeplot(dsol,[x,f3(x)],0..1);
#All 3 functions are dead zero: Not surprising with the homogeneous boundary conditions.

So the parameter you are talking enters where?
May we assume that you want a nontrivial solution?
You could experiment with giving an approximate solution as in this rather crude attempt:
dsol:= dsolve(SYS union bcs union {bcs1}, numeric,
       approxsoln=[f1(x)=100*x*(1-x),f2(x)=100*sin(Pi*x),f3(x)=100*x*(1-x)]);
plots:-odeplot(dsol,[[x,f1(x)],[x,f2(x)],[x,f3(x)]],0..1);

The particular output here I wouldn't trust right away.

With the assumptions x>0 and t>0 we can do:
restart;
sys := [v*diff(u(x,t), x) + diff(u(x,t), t) = 0, u(x,0) = exp(-x), u(0,t) = sin(t)];
res:=pdsolve(sys[1]);
eval(res,{t=0,u(x,t)=exp(-x)});
eval(res,{x=0,u(x,t)=sin(t)});
#Thus we get:
_F1:=y->piecewise(y<0,exp(v*y),sin(y));
#So that the solution is
res;
sol1:=collect(convert(res,Heaviside),Heaviside);
#sol1 is the same as
sol:=u(x, t) = exp(t*v-x)+Heaviside(t-x/v)*(sin(t-x/v)-exp(t*v-x));
simplify(sol-sol1) assuming t>x/v;
simplify(sol-sol1) assuming t<x/v;