What impresses me with example 2 is the speed. I compared with my homemade and thus amateurish DDE solver. I'm happy that it produced the same graphs as dsolve did, but as I expected it took quite a lot longer: dsolve on my computer took roughly 5 seconds to compute sol(300), whereas my program which produced a 30201x7 matrix took approximately 70 times as long.

@Rouben Rostamian  Yes, you can use an approximate solution. In this case we know the solution so the following is cheating to some degree:
##Verifying an exact solution
odetest(y(x)=1/2*sin(2*x),eval({de,bc},omega2=4));
## I'm not using the exact solution but something vaguely similar, and it happens to work:
res:=dsolve({de,bc}, numeric,approxsoln=[y(x)=sin(2.4*x),omega2=3]);

@kyanashayan If possible include the integral in its differentiated form in the system.
I will give a simple example where I use two different methods:

restart;
ode:=diff(y(x), x, x) = -(8*omega*(1-exp(-8*x))*exp(-8*x)/(2*x^2)-Pi^2/(32*x))*y(x);
bcs:=y(0)=0,D(y)(0)=1,y(1/2)=0;
## First method using 'output=listprocedure':
res:=dsolve({ode,bcs},numeric,method=bvp[middefer],output=listprocedure);
plots:-odeplot(res,[x,y(x)],0..1/2);
Y:=subs(res,y(x)); #Fishing out the procedure for computing y(x)
A:=Int(sin(x*y(x))+cos(x),x=0..1/2); #The simple example of an integral
evalf(eval(A,y=Y)); #Replacing y by Y and using numerical integration
## Second method for the same integral:
ode2:=diff(a(x),x)=sin(x*y(x))+cos(x);
res2:=dsolve({ode,ode2,bcs,a(0)=0},numeric,method=bvp[middefer]);
res2(0.5);

@Rouben Rostamian  Replacing omega^2 by omega2 solves the problem.

@farzane There is no way that you can get symbolic results in that case. You will have a characteristic polynomial of degree 8 containing 30 (!!!) parameters. As there is no general formula for the roots of a polynomial of degree 5 and above you stand absolutely no chance of getting a symbolic result.

Incidentally the exponential function exp(x) is written like that, NOT like this e^x.
You can correect that after the fact by doing as follows:

indets(Sys,`^`); #Just to see the problem
op(e^x); #The operands of e^x
Sys2:=evalindets(Sys,`^`,s->exp(op(2,s)));
indets(%,function);
Var:=[b1,c1,a2,b2,c2,d2,a3,d3];
indets(Sys2,name); #Notice the indexed v
nops(%) - 8; ## 30 (!!!)
(A, w) := LinearAlgebra[GenerateMatrix]( Sys2, Var ); #w :You cannot use v as it appears in A

op(1,A); #8x8

@al-faraj Did you try the code I gave? If you did you should have received the answer

10.11073088

And I don't think that I misunderstood your code.
I wasn't aware of a Danish version of Maple.
In your Maple version is the decimal separator (or whatever the term is) a comma?
What comes out of evalf(Pi); in your version
3.141592654
or
3,141592654

?

@maple fan The branches for a^z are those which are consequences of its definition:

a^z = exp(z*log(a)), where log is any branch of the logarithm.

Those are

log(z) = ln(abs(z)) + I*n*argument(z)

where n is any integer and where argument(z) is any argument for z.

@mskalsi Although the use of the files maple.hdb, maple.ind, maple.lib some (many) versions ago was replaced by archive files ending in .mla the use of the old method still works.
I cannot test it on the ones you mention since it seems you have to pay to get access to them.
However, I happen to have some old maple.ind, maple.lib (no maple.hdb). They contain among other things a package called DMat.

On a computer with Maple 2015 I did
libname:="E:/DMat/libDMat",libname;
with(DMat);
##The content showed up.
## In libDMat are just the files maple.ind and maple.lib

After that on another computer with Maple 12 I did the similar thing (the path is different). It worket too.

So I cannot tell you why you cannot make it work.

#Just an afterthought: Is Rif from the link really a package or is it just a collection of procedures not made into a package. If that is so then with(Rif); would result in an error.
Another thing: Is the "R" in Rif really capitalized? Should it be rif?
Thirdly, as Markiyan is mentioning, rif is implemented in Maple. In Maple 12 (I don't have access to Maple 13) when I do
?rif
up comes a help page "Overview of the Rif Command Set Version 1.1"  (just as in Markiyan's comment).

@mskalsi Are you sure that the path is correct?
My path
"E:/RIFfolder"
was just an example. You need to give the exact path to the place where you put the Rif-files.
(In particular: It doesn't seem likely that the drive is E as in my example).

@Mekai Actually I doubt that this package (Gym) is not used although I cannot be sure since the worksheet doesn't work for you either.
Anyway, what is the intended meaning of the input, which at my computer gives an error:

Error, invalid neutral operator

If I convert to 1D input I get the horribly looking

4⟦A⟧120∠+2.5⟦A⟧60∠

which surely makes no sense whatsoever (to me).

@Mekai In your very short worksheet you are trying to load a package called Gym. That package doesn't come with the standard installation of Maple. I guess it is provided by your instructor?

Do you have any reason to believe that there is a nontrivial solution?

Certainly the zero solution is a solution:

odetest({f1(x)=0,f2(x)=0,f3(x)=0,f4(x)=0},dsys3); #The trivial solution is a solution!

Maybe you should look at the link given above.

Your immediate problems can be handled by differentiating the first ode twice, the second once, keeping the remaining two as they are. The resulting system can be solved for the highest derivatives: So convertsys will work.
You have to add the extra boundary conditions you get from the first ode and its derivative, and from the second ode.

However, you will get the trivial solution.

#On reviewing this I realized I misread the highest diff order of f4: I had it at 4, it is 2.
Thus the method above can be simplified to:
Differentiate the second ode once and consider the system with that and the remaining 3 equations. Use the second ode (undifferentiated) at x=0 as an additional boundary condition. The result is still just the zero solution as expected.

This reminds me of the discussion in

http://mapleprimes.com/questions/204986-Dsolve-With-Unknown-Parameter#answer218470

So may be you ought to say where you got the ideas? It may in fact help the person, who helped with a previous situation of the same sort in finding the relevant worksheet used at the time.

You don't explain why you are doing what you are doing: e.g. why you don't just try to use dsolve on dsys3?
Are people supposed to guess that or find out themselves the hard way?

@iman The second derivative with respect to the first variable of w at the point (0,z) is written
D[1,1](w)(0,z)

Similarly with the other derivatives in IBC3. Thus IBC3 should be

IBC3:={D[2](Phi)(x, h/2)+2*D[1,1](w)(x, h/2) = 0, D[2](Phi)(x,h/2)+2*D[1,1](w)(x,-h/2) = 0, -7*D[1,1](w)(0, z)+2*Phi(0, h/2)-Phi(0, -h/2) = 0, -7*D[1,1](w)(L, z)+2*Phi(L, h/2)-Phi(L, -h/2) = 0, Phi(0, z) = 0, Phi(L, z) = 0, u(0, z) = 0, u(L, z) = 0, w(0, z) = 0, w(L, z) = 0};

Now the next problem you face is not syntax, but the fact that the algorithm used cannot handle derivatives that are not normal to the boundary: The error message says as much:

In IBC3 the derivative D[2](Phi)(x, h/2) is perfectly fine: It is a derivative w.r.t. to the second variable (z) and the boundary is z=h/2, thus that derivative is normal to the boundary.
However, the derivative D[1,1](w)(x, h/2) is not normal to the boundary as it is a (second) derivative w.r.t. the first variable (x).
Maybe you intended that derivative to be D[2,2](w)(x, h/2) ?
But your problems don't end there: IBC3 also has Phi(0, h/2), which is just Phi at one point. You need to give Phi along an axis, i.e. the different ibcs should each depend on exactly one of the two independent variables (x or z).

Finally, if all that is taken care of, you will most likely run into a problem with the appearance of diff(Phi(x, h/2), x, x) in your system of pdes.
As a simple example of the problem:
restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x,x); #The heat equation
res:=pdsolve(pde,{u(x,0)=x,u(0,t)=0,u(1,t)=0},numeric);
res:-animate(t=1); #Perfectly fine
#Now consider the following changed version:
pde2:=diff(u(x,t),t)=diff(u(x,t),x,x)+diff(u(x,0),x); #Notice the appearance of diff(u(x,0),x)
res:=pdsolve(pde2,{u(x,0)=x,u(0,t)=0,u(1,t)=0},numeric);
res:-animate(t=1); #ERROR

@golnaz I think that I actually answered that in giving

plot([2*IB(z,1/5,1/2),z,z=0..1]);

You just have to replace IB by
IBcrt:=(t,a,b)->IB(t^(1/3),a,b);

as I actually implied above in giving IBsqrt.
I have noticed that you changed t^(0.5) fo t^(1/3), so I changed IBsqrt to IBcrt.