I don't see any attachment.

I just saved the attached worksheet to the cloud first as public (but that takes a while to show up and probably won't as it is not of general interest) and then as private.
The private one was retrieved and looked like the attached worksheet, i.e. in good shape.
MaplePrimes17-08-23_vokaltest.mw

I realize that I didn't use the workbook format. My test was with a simple worksheet.

So you have a concrete ode of order three, which Maple could easily solve using dsolve:
 

restart;
ode:=diff(y(x),x,x,x)+2*diff(y(x),x,x)-9*diff(y(x),x)-18*y(x)=18*x^2-18*x+22;
dsolve(ode);

Now what is f supposed to represent?
I ask because if you want to write this ode as a first order system (a very common approach), you need to introduce 3 functions f, g, and h (or whatever you would like to call them) with f(x)= y(x), g(x) = diff(y(x),x), and h(x) = diff(y(x),x,x).

Besides, you are talking about an ivp (i.e. an initial value problem), but what are the initial values?
## Solving with initial values using ics as given here:

ics:=y(0)=y0,D(y)(0)=y1,(D@@2)(y)(0)=y2;
dsolve({ode,ics});

I don't think that the idea is viable. Patches are best left to professionals.
We happy amateurs can once in a while find a workaround for some bug or inadequacy. In fact that can be a lot of fun.

Do you envision a version number on, say, patch12? Obviously, the patch idea if successful would make patch12 grow and need modification all the time. By whom? Anybody?
Is your vision a Maple 12 with no errors although with no improvements or features from later versions?

Remember that the update of the Physics package wasn't done by an amateur!
 

You have missing multiplication signs in several places in eq1, eq2, eq3.
When I execute eq1 as you wrote it I get (copying as an image):

If I insert multiplication signs where they are clearly needed I get the input equations:
 

eq1:=(diff(psi(x), x))^2+(diff(u(x), x)+4*(diff(w(x), x))^2)*(diff(psi(x), x))^2+3*(diff(diff(w(x), x), x))+5*(diff(diff(w(x), x), x))*(diff(psi(x), x))-7*(diff(diff(diff(u(x), x), x), x))-28*(diff(diff(w(x), x), x))^2-(21/2)*(diff(diff(diff(w(x), x), x), x))*(diff(w(x), x))+3 = p;
##
eq2:=20*(diff(diff(psi(x), x), x))+50*(diff(diff(diff(w(x), x), x), x))+8*(diff(diff(diff(diff(psi(x), x), x), x), x))-7*(diff(w(x), x))+7*psi(x) = 0;
##
eq3:=4*(diff(diff(w(x), x), x))-4*(diff(psi(x), x))+34*(diff(diff(diff(psi(x), x), x), x))+26*(diff(diff(diff(diff(w(x), x), x), x), x)) = 0;

As you will notice, they are quite different from yours. That is a consequence of the bad decision (in my opinion) of allowing spaces to be interpreted as multiplication in 2D math input. People forget the space occasionally.
I never use 2D math input.
I shall correct my answer below to reflect the revised equations.
I didn't see the problem till today. And why not? Because the original equations also make sense and happen not to become totally ridiculous (that would have saved us).
Example:
In the original eq3 you have a missing multiplication sign in the term (23+11)(diff(psi(x), x, x, x)).
Now this is automatically simplified to 34(diff(psi(x), x, x, x)). In Maple numerical constants like (here) 34 act as constant functions so that 34(y) = 34 for all y no matter what y is. Thus 34(diff(psi(x), x, x, x)) = 34.

@torabi I chose the condition (D@@2)(u)(0) = 0 as the first simple one I could think about.
The quantities that can be involved in the boundary conditions for your system are all derivatives of orders from j = zero to one less than the highest appearing for each function.
Thus the quanties that can be involved in boundary conditions for your system are:
 

N:=[4,3,4]:
nms:=[psi,u,w];
pts:=[0,1];
S:={seq(seq(seq((D@@i)(nms[j])(a),i=0..N[j]-1),j=1..nops(N)),a=pts)};

Now you already have 11 boundary conditions and they are all of the simple form (D@@i)(f)(a) = c for some i ,f, and c.
So if you want to add an extra simple condition of the same form you can do:
 

## The bcs you have given we shall name bcs:
bcs:={psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0};
## So you can choose from
S minus lhs~(bcs);

Pick one of those and equate it to some concrete value c, add it to bcs, and use dsolve.
Be aware that p will depend on what you end up doing!
########## To illustrate the last remark try all of this:
 

S1:=S minus lhs~(bcs);
for s1 in S1 do
  try
    res[s1]:=dsolve({eq1,eq2,eq3} union bcs union {s1=0},numeric)
  catch:
    printf(cat(StringTools:-FormatMessage(lastexception[2..-1]),"\n"))
  end try
end do;  
indices(res); # The successes
## We have revised the p results to reflect the revised eqs:
## In all 3 cases we just get the trivial solution
res[(D@@2)(u)(0)](0); # p = 3
res[(D@@2)(u)(1)](0); # p = 3
res[D(u)(0)](0); # p = 3
## Now set the quantities in S1 equal to 1 instead of 0:
for s1 in S1 do
  try
    res1[s1]:=dsolve({eq1,eq2,eq3} union bcs union {s1=1},numeric)
  catch:
    printf(cat(StringTools:-FormatMessage(lastexception[2..-1]),"\n"))
  end try
end do;
indices(res1); # The successes (the same as before)
res1[(D@@2)(u)(0)](0); # p = 13.5000000390248 (the revised eqs)
res1[(D@@2)(u)(1)](0); # p = -17.9999999849048 (rev. eqs)
res1[D(u)(0)](0); # p = -39.0000000347809 (rev. eqs)

If you do the same loop with s1 = 2 you get the following p-values: p = 24, p = -39, p = -81.
So you should ask yourself what you are trying to do!
## If you want to see the graphs of e.g. u in the different cases you can do:

plots:-display(Array([seq(plots:-odeplot(res[rr],[x,u(x)]),rr=indices(res,nolist))]));
plots:-display(Array([seq(plots:-odeplot(res1[rr],[x,u(x)]),rr=indices(res1,nolist))]));

@torabi I don't know what you mean by " is it possible for solving boundary condition.? ".
You haven't commented on my statement " If your intention is to determine p at the same time you solve ...."
I take that to mean that you do want to determine p from the boundary conditions. But then you must come up with an additional boundary condition.
As an example I added the condition (D@@2)(u)(0)=0 besides the one you just gave above: D(u)(1)=0:
 

dsys3a := {eq2, eq3, eq1, psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0,(D@@2)(u)(0)=0};
res:=dsolve(dsys3a,numeric);
res(0); # you will see that  p = 3.00000000000000. #Using the revised equations: See above
plots:-odeplot(res,[x,u(x)],0..1); #Notice that u is very, very small.

The solution for u, psi, and w above is really the trivial one u = w= psi = 0 identically.
If on the other hand you don't want to impose other boundary conditions than the 11 you got, then you must supply the value of p as I have done in these lines where  I took p = 1:
 

dsys3b := {eq2, eq3, eval(eq1, p = 1), psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0};
res1:=dsolve(dsys3b,numeric);
plots:-odeplot(res1,[x,u(x)],0..1);

@Christopher2222 You could use the following script in your maple.ini file:
 

ModGasDensity:=proc(elem) local vuu;
  uses ScientificConstants;
  if not HasElement(elem) then error "Element %1 is not included", elem end if;
  vuu:=rhs(GetElement(elem,'density')[2]);
  if convert(subs(vuu,units),`global`) = 'kg/L' then
    vuu:=[vuu[1],vuu[2],'units'='g/L']; 
    printf("Modifying the unit for the density of the gas %a from kg/L to g/L\n",elem);
    ModifyElement(elem,'density'=vuu)
  else
    NULL
  end if
end proc:
for i from 1 to 110 do 
   try
     ModGasDensity(i)
   catch:
   end try
end do;

You may want to comment out the printing.
Test:

with(ScientificConstants);
GetValue(Element(Krypton,density(gas)));
Units:-UsingSystem();

### Since Maple 8 through 15 have L where Maple 16 and up have kg/L I modified the script some and included it in a procedure call to avoid assignments to global variables.
This should work in all versions (I hope). At least it works in Maple 12, 15, 16, 2017.
In Maple 8 we cannot use 'uses ScientificConstants' so the long form have to be used for HasElement and the rest. I have done that, but another problem with 8 is that it doesn't read the maple.ini file in my Users folder. Not a big problem since I only look at Maple 8 for purposes of comparison with later versions.

proc() local VERSION,ns,N, ModGasDensity,i;
   VERSION:=convert(kernelopts(version),string):
   ns:=StringTools:-Search(",",VERSION):
   N:=parse(VERSION[7..ns-1]);
   if N<8 then return NULL end if;
   ModGasDensity:=proc(elem) local vuu,U;
     if not ScientificConstants:-HasElement(elem) then error "Element %1 is not included", elem end if;
     vuu:=rhs(ScientificConstants:-GetElement(elem,'density')[2]);
     U:=convert(subs(vuu,units),`global`);
     if N>=16 and U = 'kg/L' or N<16 and U='L' then
       vuu:=[vuu[1],vuu[2],'units'='g/L'];
       printf("Modifying the unit for the density of the gas %a to g/L\n",elem);
       ScientificConstants:-ModifyElement(elem,'density'=vuu)
     else
       NULL
     end if
   end proc;
   for i from 1 to 120 do
     try
       ModGasDensity(i)
     catch:
     end try
   end do;
end proc();

@Mac Dude Right, it ought to be easy to fix.
To see the units used in the record it is easier to do
ScientificConstants:-GetElement(H,density);
    

In Maple 8, 12 and 15 (probably 9, 9.5, 10, 11, 13, and 14 too) the last item is units = L.
So yes, certainly the part of ScientificConstants having to do with the density of gases couldn't have been used much over the years: surviving from Maple 8 (when the package was introduced) to Maple 15 with a density reported as 1000 times too low, and then from Maple 16 to Maple 2017 with a density that is 1000 times too large (both when SI units are desired).
 

A comment about the unit used at the time the data were put into the system:
In Maple 15 the unit reported by line 42 in the procedure GetValue  by
ScientificConstants:-numTabElem[op(x)[1]]:-props[op(0,op(x)[2])]:-units
is L, i.e. liter, a measure of volume! The 0.088L is then correctly converted to the SI unit m^3 as 0.000088.
##
In Maple 16 (and following versions) this should have been corrected to the unit g/L, but was made kg/L, resulting in 88.0.

## Note x g/L = x kg/m^3.

@9009134 You have a problem with a wide range in the numerical coefficients:
Take a look at
 

(I set Digits:=15)
sys3;
SYS;
## The absolute values of the floats in the systems:
abs~(indets(evalf(sys3),float));
(min,max)(%); # Answer: 4.50605975510204*10^(-80), 4.66666666666667*10^7
abs~(indets(evalf(SYS),float));
(min,max)(%); # Answer 1.32275132275132*10^(-41), 4.82434787943840*10^99

I don't know if you will be able to scale this problem so that this huge relative difference in size disappears or at least gets considerably smaller. Otherwise you will have to work with a large value of Digits, way into the software floating point range, i.e. much larger than 15, which will slow down the handling considerably.

@9009134 You must make sure that the different quantities actually get defined before embarking on the dsolve process.
When I execute your worksheet ss.mw by using !!! in the top panel I find that SYS evaluates to NULL, i.e. the solve command has no solution. Furthermore newsys2 evaluates to

which is not a system of odes.
So my suggestion is to write the whole thing in a worksheet (not a document) using 1D math input (aka Maple input). 
Write one command at a time and examine the output to see if there is anything wrong at each step.
As it is right now it is a mess (quite frankly).

@Dmitrii It is highly unlikely that you would want a function just defined on your list x having corresponding values in your list y. But just in case you do, that is easily done:
 

restart;
x:=[0,1,2,3,4,5];
y:=[0,2,4,6,8,10];
for i from 1 to numelems(x) do f(x[i]):=y[i] end do;
eval(f); #f has option remember
op(4,eval(f)); # The remember table
f(t); # returns unevaluated.
f(0.5); # returns unevaluated since f(0.5) was never defined.
f(2); # returns 4 since 2 is in the remember table.

Surely, you want some kind of interpolation as Rouben is saying.

If you are pretty sure that this behavior applies to any Maple product I would very much like to see the example in Maple itself.
Give us text or upload a worksheet, please.

@John Fredsted Maybe Tom Leslie's example was too simple to show the problem clearly.
Anyway, I made a slightly more complicated example, where the first method of solving successively is fine, but the second is not and should be avoided:
 

restart;
ode1:=diff(x(t),t)=x(t);
res1:=dsolve(ode1);
ode2:=diff(y(t),t)=y(t)+x(t);
res2:=dsolve(eval(ode2,res1)); #Fine
## And is the same as solved as a system:
dsolve({ode1,ode2}); 
########################
## This way is no good:
dsolve(ode2,y(t)); 
eval(%,res1);
value(%);

@tomleslie I wouldn't complain about your first example either. It amounts to
 

restart;
de1:=diff(y(x),x)=1;
de2:=diff(u(x),x)=1
dsolve(de1);
dsolve(de2);

In both cases the arbitrary constant is called _C1. That is well known and can be counted on.
Consistency is important. In other parts of Maple constants are renumbered every time you redo the calculation as in _Z3 and next time _Z4 etc. which is rather annoying.

to 10 do solve(sin(x)=0,allsolutions) end do;