I have almost no experience with the discrete fourier transform, but maybe that is what you need:

restart;
eq1:= diff(l(t),t) = Pl(t)/(0.22124+0.2);eq2:= diff(Th(t),t) = PTh(t)/ (0.2 * l(t)^2);
eq3:= diff(Pl(t),t) = PTh(t)^2/(0.2 * l(t)^3) - 0.22124*9.8 + 0.2*9.8*cos(Th(t));
eq4:= diff(PTh(t), t) = - 0.2 * 9.8 * l(t) * sin(Th(t));
eq5:= F(t)= 0.22124 * 9.8 + 0.2 * diff(Pl(t),t)/(.42124);
ans:= dsolve({eq1, eq2, eq3, eq4, eq5, l(0)=0.449, Th(0)=.94, Pl(0)=0, PTh(0)=0}, [l(t), Th(t), Pl(t), PTh(t), F(t)], numeric,output=listprocedure);
#I shall work with the fourier transform of F.
Fn:=subs(ans,F(t));
#Aping the help page for FourierTransform:
with(DiscreteTransforms):
Digits:=15:
Z:=Vector(1001,i->Fn(-5+(i-1)/100),datatype=complex[8]);
Z2:=FourierTransform(Z);
#Checking by going back:
Z1:=InverseFourierTransform(Z2);
#We get back to Z OK:
fnormal~(Z-Z1)[1..10];
#Corresponding times:
T:=Vector(1001,i->(-5+(i-1)/100),datatype=float[8]);
#simplify removes the ...+0.*I terms in Z:
Z[1..10];
p0:=plot(T,simplify(Z),color=blue): p0;
p1:=plot(T,simplify(fnormal~(Z1)),color=red): p1;
plots:-display(p0,p1);
#About plotting the discrete fourier transform Z2 I don't know what is expected. Z2 is a vector with complex entries.

#####################
A comment:

By uniqueness this shows that Pl and PTh are odd,  and l, Th, and F are even (initial conditions also taken into account):
PDEtools:-dchange({Pl(t)=-Pl2(tau),l(t)=l2(tau),PTh(t)=-PTh2(tau),Th(t)=Th2(tau),F(t)=F2(tau),t=-tau},[eq1,eq2,eq3,eq4,eq5],
[tau,Pl2,l2,PTh2,Th2,F2]);
subs([tau=t,Pl2=Pl,l2=l,PTh2=PTh,Th2=Th,F2=F],%);
#Compare with:
[eq1,eq2,eq3,eq4,eq5];

You are doing nothing wrong. Maple needs some help.

E:=eliminate({eqn1, eqn2, eqn3},{H,Tf});
#A start guess for fsolve can be found by plotting:
plot(E[2],Tr=40..100);
fsolve(E[2],Tr=45); #It helps with a start value
#Interestingly enough also this works
fsolve(E[2],Tr,complex);
eval(E,%);

## In fact fsolve can find this solution at once like this:

fsolve({eqn1, eqn2, eqn3}, {H, Tf, Tr},complex);

I don't know why you reposted the problem, see

http://www.mapleprimes.com/questions/145690-How-To-Set-The-Boundary-Conditions#comment145700

but I still think that the quote I gave from the help page for pdsolve/numeric in my comment to your previous question tells the story: It is a time based solver. If that means what I think it means, then it starts at some time t0 with given values and works its way forward in time and on the way also obeying whatever boundary conditions might have been given. You don't have any initial conditions for it to start with.
Notice that I'm not saying anything about the existence or uniqueness of a solution to your problem. I'm just pointing out the limitations of the time based solver.

Added: You have a condition ui(0, t) = 0. Did you intend that to be ui(e,0) = 0? If so, then you have one initial condition, but you need one for ur too. And then you have too many boundary conditions.
Furthermore, notice that the zero solution ui(e,t) = 0, ur(e,t) = 0 satisfies your pdes and all your conditions.

In your worksheet http://www.mapleprimes.com/ViewTemp.ashx?f=0_1365331923/code_for_fdm.mw

you could insert the line

V:=Array(-1..M+1,0..N);
before assigning to V in the loops.

In your original version V is a table (although implicitly created).

Arrays are easier to deal with.

However, be aware that unassigned elements will be zero. So before doing it try executing the table version (the original) and then do these lines:
indets({entries(V,nolist)},name);
Ve:=map(eval,V):
indets({entries(Ve,nolist)},name);
Notice the change from V to Ve.
But since also Ve has unassigned V-elements you should examine your loops.

Since  x<0 and x>1 is never satisfied, the 'otherwise' takes effect (by default 0).
Try specifying a default (e.g. 2). Here just plotting:

plot(piecewise(0 > x and x > 1, 0, x > 1, 1, 2), x=-5..5);

What did you actually want that piecewise to look like?

The student[intersect] was not meant for doing the job, which (in Maple 17) algcurves[intersectcurves] is doing and what is (apparently) what you want.

In Maple 17

algcurves[intersectcurves](C1,C2,x,y);

results in

[[7, [x, y, 1]], [2, [1, 0, 0]]]

Since obtaining the answer involves finding the roots of a polynomial of degree 5 you are almost surely not able to get any better result than a result involving RootOf. 
RootOf(pol) stands for any of the roots of the polynomial pol (here in the variable _Z introduced by isolate).

indets(M1sol,RootOf);
pol:=op([1,1],%);
#You can use an alias for RootOf(pol):
alias(alpha=RootOf(pol));
M1sol;

##I add a concrete example:

nm:=indets(M1alone,name) minus {M1};
ex1:=eval(M1alone,nm=~{i$i=1..nops(nm)});
isolate(ex1,M1);
evalf(%);

#You will notice an index in the RootOf.
Since sqrt is a specific square root (the principal) it is to be expected that not all the roots of the polynomial will give rise to a solution of the original equation.

Well, it says in ?algcurves,intersectcurves

that

"The algcurves[intersectcurves] command was introduced in Maple 17."

which implies that it is absent from earlier versions.

restart;
Normalizer:=radnormal;
assume(s>0,s<1);
r:=2*arctanh((1+s)*tanh(x)/sqrt((1-s)*(1+s)))+ln((1/2)*arccosh(1/s)-x);

L:=limit(r,x=1/2*arccosh(1/s));
simplify(L);
radnormal(%,rationalized);
combine(%);
radnormal(%);

I have allowed myself the assumptions s> 0, s < 1, beta>0:
restart;
Normalizer:=radnormal; #Makes the result look better.
r1:=2*(1-s^2)^(1/2)*arctanh((1+s)*tanh(x)/((1-s)*(1+s))^(1/2))/((1-s)*(1+s))^(1/2)+ln((1/2)*arccosh(1/s)-x);
x:=1/2*arccosh(1/s)-beta*u;
r1;
series(r1,u,2) assuming s>0,s<1,beta>0;
res:=simplify(%) assuming s>0,s<1,beta>0;
#Numerical evidence for the correctness:
plot(eval([convert(res,polynom),r1],{beta=1,s=1/2}),u=0..0.5);

#I have edited my answer to reflect the change made in the comments.

Try laplace transformation. You can use that laplace works easily with convolutions. Your integral is a convolution of f with a 3rd degree polynomial.
So in Maple:
restart;
with(inttrans):
u:= f(t) = exp(-3*t)+t^3+Int((7-x+(17/2)*x^2-(35/3)*x^3)*f(t-x), x = 0 .. t);
laplace(u,t,s);
isolate(%=0,laplace(f(t),t,s));
invlaplace(%,s,t) assuming t>0;
ft:=rhs(%):
plot(ft,t=0..1);

@dey_meghdeep The computation of the eigenvalues and eigenvectors (as you present it) involves using all the exact roots of the characteristic polynomial of dyna. These roots can only be presented in the form of RootOf with indices 1 .. 21. The output for L is huge (without evalf).
Using evalf on Eigenvectors(dyna) with Digits = 10 as you are doing is likely to create a substantial accumulation of roundoff errors.
Instead of using evalf with a higher setting of Digits (e.g. evalf[15](Eigenvectors(dyna)) I would suggest letting dyna be a floating point matrix. This could be done by changing dyna:
dyna:=Matrix(21,dyna,datatype=float[8],shape=symmetric);
I have given it shape=symmetric (well it is) to avoid getting numerical results of the type 1.234+0.*I, although these could easily be removed by simplify.
If you want to see the results use
interface(rtablesize=infinity);  #or just interface(rtablesize=21);

It is hard for me to see the images. But as far as I can see your summation over j goes from 1 to 3, so it is a sum of a concrete number of terms.
When you say that you want to differentiate w.r.t. P[j,k] and you are summing over j, I cannot follow you.

An example (since I cannot see clearly what you got):

S:=sum(f(P[j,k]),j=1..3);
#It makes sense to differentiate w.r.t. any particular P e.g. P[2,k]:
diff(S,P[2,k]);
#And even w.r.t P[j,k] with j unassigned, but then of course the answer will be 0.

There is no problem with the first one as you noticed yourself (why didn't you give us Maple code instead of Mathematica code and in text form?).
For the second model use options maxfun (maximal number of function evaluations) and abserr (absolute error).
See ?dsolve,numeric and ?maxfun.

restart;
sys:={m*diff(x(t),t,t)=spring(x(t))+friction(diff(x(t),t)),x(0)=1,D(x)(0)=0};
m:=1:
beltv:=.1:
spring:=x->1000*(1-x);
viscous:=v->-30*(v-beltv);
friction:=viscous;
res1:=dsolve(sys,numeric);
plots:-odeplot(res1,[t,x(t)],0..1);
plots:-odeplot(res1,[[t,diff(x(t),t)],[t,beltv]],0..1);
#Second system:
coulomb:=v->-25*signum(v-beltv);
plot(coulomb,0..0.2); #The discontinuity of coulomb is the problem
friction:=viscous+coulomb;
res2:=dsolve(sys,numeric,maxfun=10^7,abserr=1e-5);
plots:-odeplot(res2,[t,x(t)],0..1);
plots:-odeplot(res2,[[t,diff(x(t),t)],[t,beltv]],0..1);

The easy way is to assign to x0, y0, and r0 instead of x[0],  y[0], r[0]. Then you will have no problem.
In Maple 17 you can use x__0 (double underscore). It prints like x[0].