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' ?

@Konstantin@ A common (and recommended) method is to include the integrals in the system.

Suppose as a simple example that you want to find the integral int(z(x)*q(x),x=0..a), where q is a known function.
Then you could include the ode diff(w(x),x)=z(x)*q(x) with w(0)=0 in the system:

sys:={diff(z(x), x$2)+z(x) = -2, diff(w(x),x)=z(x)*q(x), z(0)=0, z(1)=1, w(0)=0}.

Then the matrix returned with the option=Array(...) will include w-values, and of course w(a) will be your integral.

Did you actually try this:

dsolve({eq,y(0)=1,D(y)(0)=0},y(x),type='series',order=10);

It works for me in Maple 18, Maple 15, and Maple 12. I'm pretty confident that it works in all versions from 12 and up. Maybe it also works even in some earlier releases. I don't have access to those on this machine.

As for series:

series(cos(x),x=0,10);