@emmantop But bd represents the fourth derivative, not the second derivative, so it cannot be used in eq[12].

To illustrate that we can do:
restart;
bd := (f[i]-4*f[i-1]+6*f[i-2]-4*f[i-3]+f[i-4])/h^4;
bd0:=eval(bd,i=2);
for i from -2 to 2 do f[i]:=(i*h)^4 end do; # Corresponding to f(x)=x^4
bd0; #Gives the correct fourth derivative of x^4.
for i from -2 to 2 do f[i]:=(i*h)^2 end do; # Corresponding to f(x)=x^2
(f[1]-2*f[0]+f[-1])/h^2; #This gives the correct second derivative of x^2
bd0*h^2; # This does not

@smith_alpha You could add the .txt by using  cat("D:/",file_name,".txt").

If you really defined f as a procedure, then plot(f,-2..2) should work.

What was your definition of f (exactly)?

@John Starrett With Digits:=30 as I suggested earlier I cannot see any difference between 
display(Array([PW1a,PW2a,PW3a])) and display(Array([PW1,PW2,PW3])).

Before computing DEP1 and DEP1a I set Digits back to 10 and used maxfun=0 in both.
Again no difference.

Note:  I used the first version of your posted worksheet.

@emmantop I would linearize after discretization. By linearization I simply mean using Newton's method on the system of equations ensuing from discretization. That involves a starting guess and iteration like in the one dimensional Newton's method.
However, I would begin by writing the ode-system as a first order system. Then discretize.

Obviously you would have to replace infinity by some number like 10 and use a midpoint method.
To test your results you could compare with

res:=dsolve({ode1,ode2} union subs(infinity=10,{bc1,bc2}),numeric,method=bvp[midrich]);

@adel-00 
The straightforward approach actually works, you just have to interpret the result correctly (in the first version of this comment I didn't have it straight).
eqs:=subs(x(t)=x,y(t)=y,z(t)=z,rhs~(dsys))=~0;
res:=solve(eqs,{x,y,z});
#x=RootOf(conjugate(_Z)+_Z) means that x is a number whose real part is 0.
#Thus x = I*x2, where x2 is real.
#The result is therefore
evalc(subs(RootOf(conjugate(_Z)+_Z)=I*x2,res));
#i.e. there are infinitely many solutions parametrized by x2.
#We check these solutions below.
#We could split instead:
eqs2:=subs(x(t)=x1+I*x2,y(t)=y1+I*y2,z(t)=z,rhs~(dsys))=~0;
eqs2R:=(evalc@Re)~(eqs2);
eqs2I:=(evalc@Im)~(eqs2);
solve(eqs2R union eqs2I,{x1,x2,y1,y2,z});
#Notice that the solutions are the same as found directly above.
#Checking:
eval(eqs2,%);
simplify(%) assuming real;
 

Please note that I revised my second answer totally. I submitted a bug report expressing this conflict about plotting before saving.

@adel-00 I think you meant 0.25*z(t)^2+Re(y(t))^2+Im(y(t))^2 because that quantity should be and appears to be (roughly) constant (=1/16=0.0625), whereas the one you actually have doesn't:

plots:-odeplot(res,[t,0.25*z(t)^2+Re(y(t))^2+Im(y(t))^2 ],0..10);
plots:-odeplot(res,[t,0.25*z(t)^2+Re(y(t))^2+Im(y(t)) ],0..10);

By saying "I have written a demonstration worksheet to show this problem", I presume that you intended to include that worksheet in your question, but I don't see any link to it.

@adel-00 If you just want to plot the quantity 0.25*z(t)^2+Re(y(t))^2 then you can just use res from before:

plots:-odeplot(res,[t,0.25*z(t)^2+Re(y(t))^2],0..2);
# but if you want to get actual numbers you can use output=listprocedure, which would be OK to use in any case.
#Notice that this time I have closed the assignments to resL, Y, Z by a colon. If you replace colon with semicolon you will know why (a huge output: it is usually hidden. It is surely not intended. I shall report it as a minor bug).

resL:=dsolve(dsys union {x(0)=0,y(0)=0,z(0)=-1/2},numeric,output=listprocedure):
#You can use odeplot as before:
plots:-odeplot(resL,[t,0.25*z(t)^2+Re(y(t))^2],0..2);
Y,Z:=op(subs(resL,[y(t),z(t)])):
Y(0.12345);
Z(0.12345);
#You could if you wish define the quantity as W:
W:=t->0.25*Z(t)^2+Re(Y(t))^2;
W(.987654);
#And you can plot W (same plot as above):
plot(W(t),t=0..2);

@RAfossey Here is a stepwise solution by laplace transform. As you are doing I'm assigning to a,b,c,d immediately although that doesn't seem to be what the problem formulation called for.

restart;

with(inttrans):
a:=1;b:=1;c:=1;d:=1;
sys:=diff(x(t),t,t)=a*x(t)+b*y(t),diff(y(t),t,t)=c*x(t)+d*y(t);
ics:=x(0)=1,D(x)(0)=0,y(0)=0,D(y)(0)=1;
#Applying laplace elementwise (notice the tilde ~ ) to the system:
laplace~([sys],t,s);
eval(%,{ics});
#Solving for the laplace transforms:
L:=solve(%,{laplace(x(t),t,s),laplace(y(t),t,s)});
#Elementwise inverse transformation
sol:=invlaplace~(L,s,t);
#Extracting the solutions:
X,Y:=op(subs(sol,[x(t),y(t)]));
#Plots:
plot([X,Y],t=0..2);
#This is a parametric plot:
plot([X,Y,t=0..2]);

@wenny You should try executing the worksheet again. I did not get that error message.
If you get that error once again then try this:

indets({g1,g2,g3},name);

you should be seeing the unassigned names and hopefully what you see are just the variables.

@ I'm apalled to learn that some math professor out there cares a hoot about how you indicate the end of proof as long as you do it!

@jonlg Maple simply cannot do the integral even when T1, T2, U0, L0 are given concrete values.

I dropped the case t<=0 because I thought (perhaps wrongly) that you had no interest in t<0.

piecewise works like this: if the first condition is satisfied then the following expression is returned. If not then the next condition is evaluated, etc. If none of these is satisfied then 0 is returned. That 'otherwise' case is specified explicitly in my example as the last argument; when given explicitly it can be anything.

Could you give us a (simple) example of what you have in mind, thereby answering questions like 'Does the parameter have to be determined also or can it chosen arbitrarily' ?