It looks to me like you're plotting the acceleration, which according to your differential equation is

And this will indeed be discontinuous, because the signum function is discontinuous.  At the local minima and maxima of x(t), where the velocity changes sign, its signum changes from +1 to -1 or vice versa.

Note that in Maple π is Pi and √-1 is I.  You might also find the mul, expand and evalf commands useful (maybe also evala).

Your "hint" makes no sense as it stands. 
All the Taylor coefficients of your f(x) about 0 are 0, so the remainder is just f(x) itself.
The remainder, according to Taylor's theorem with the Lagrange form of the remainder,
should be (x-a)^(n+1)/(n+1)!*(D@@(n+1))(f)(xi) for some xi between a and x.  My guess is that
you are supposed to evaluate this for several different n's, with xi = 0.25 and x = 0.5,  and plot the points
[n,  (0.5-0)^(n+1)/(n+1)!*(D@@(n+1))(f)(0.25)] you obtain.  However, there won't be a limit as n -> infinity:
the xi's that make (x-a)^(n+1)/(n+1)!*(D@@(n+1))(f)(xi) equal to the remainder will be different for different n.

It looks to me like the problem is a bug in 2-D math with restart and alias.

In a Standard worksheet, using 2-d math input:

> restart;
   alias(J=BesselJ);

Now click the !!! button repeatedly.  On the odd-numbered clicks the result is J (as it should be), but on the even-numbered clicks it is nothing (and J has not been aliased to BesselJ).

The problem does not occur if you use 1-d Maple input, or if you put the alias command in a separate input group.

Yes, of course it has numerical solutions.  For any x you can use fsolve to find them.  You can also plot the curve, either using plots[implicitplot] or algcurves[plot_real_curve].  For the latter, note that the equation (at least for y > 0) is equivalent to (y^2 - (1+x)^3 - (1+x)^4)^5 = y.

> eq:= (y^2 - (1+x)^3 - (1+x)^4)^5 - y;
   algcurves[plot_real_curve](eq, x, y);

Sometimes the crashes seem random, but this one seems to consistently crash Maple
(both Maple 13 and Maple 14) on Sun SPARC Solaris.

>    with(Statistics):
      n:= 19:
      X:= <seq(sin(1.0*j),j=1.. n)>;
      w:= <seq(j, j=1..n)>;
      Median(X, weights=w);

The crash occurs in the external procedure MapleMedian called by Statistics:-DescriptiveStatistics:-Median:-GetValue.

Under Windows 7 that sequence of commands did not usually crash Maple 14 Classic, but the following one usually did:

restart; with(Statistics):
       for n from 17 to 19 do
       X:= <seq(sin(1.0*j),j=1.. n)>;
       w:= <seq(j, j=1..n)>;
       print(n,Median(X, weights=w));
 end do:

print is a function.  Its argument needs to be in parentheses.  So

for i from 1 to 200 do
  if ithprime(i)=prevprime(500) then print(i); end if;
end do;

should work fine. 

The attempt using solve is doomed because Maple doesn't know how to solve equations involving ithprime (and it shouldn't, because ithprime is only defined on positive integers).

But for a more efficient way, see numtheory[pi].

You have 6 equations (because of symmetry) and 9 unknowns.   Thus there are 3 extra degrees of freedom: in principle you should have solutions where 3 of the unknowns are arbitrary, the rest expressed in terms of them.
Writing Sij for Psi[i](j), your equations will be of the form

> eqs:= [S21*S03 - a1, S01^2+S11^2+S21^2-2, S02^2+S12^2+S22^2-2, S03^2+S13^2+S23^2-2,
 S11*S02 + S21*S12 - a2, S12*S03 + S22*S13 - a3];

where a1, a2, a3 are functions of T.  You could then try e.g.

> solve(eqs, {S01, S02, S03, S11, S12, S13, S20, S21, S23});

I don't know which of the solutions returned will produce real values.

The basic problem, I think, is that the curves F(x,cp) = 0 come very close to the curves where F(x,cp) -> +/- infinity, so that some of the regions F(x,cp) < 0 or F(x,cp) > 0 are very narrow.  Unless you happen to get grid points falling in those narrow regions, you miss part of the curve.  This can be seen more clearly in

> plots[implicitplot](tan(x)-y, x=0..2,y=-200..200, gridrefine=3,signchange=false,grid=[10,10]);

Note that even though I asked for y = -200 to 200, only the part from about y=-80 to 65 is shown.  For y = 65, the region where tan(x) - y > 0 is only from x = arctan(65) = 1.555412925 to Pi/2 = 1.570796327, and it gets even narrower as y increases.

Sorry, Maple has given you your answer.  The current facilities for numerical PDE solving only allow two independent variables.  Unless there's a way to reduce your problem to something easier to solve, you'd have to write your own PDE solver. 

Assuming some of those products really are nonzero (and you're using Multiply in the LinearAlgebra package rather than "multiply"), that ought to work.  Maybe if you uploaded your worksheet we could spot the problem.

The problem is your boundary condition on the second derivative.  Try (D@@2)(f)(2) instead of
(diff(f(eta), `$`(eta, 2)))(2)

Maple has trouble with these Bessel integrals.  For example,

int(exp(k*cos(t))*cos(n*t),t=0..Pi) assuming k>0, n::integer;

works (i.e. produces Pi*BesselI(n,k)), while if you replace the Pi by 2*Pi or the 0 by -Pi the integral returns unevaluated.   So what you can do is

1) use IntegrationTools[Change] to rescale x so that the integral is from 0 to 2*Pi.

2) use evalc to separate real and imaginary parts.  Note that by symmetry the imaginary part is 0, while the
real part is 2 times an integral from 0 to Pi that Maple can do.

Not directly.  Expand your f(theta) in a Fourier cosine series

f(theta) = c[0]/2 + sum(c[j]*cos(j*(theta-theta0)*Pi/beta), j=1..infinity)
and the solution is (if I haven't made an error)

A(r,theta) = c[0]/2 + sum(c[j]*((r/R1)^(j*Pi/beta) + (R1/r)^(j*Pi/beta))
/((R2/R1)^(j*Pi/beta) + (R1/R2)^(j*Pi/beta)) * cos(j*(theta-theta0)*Pi/beta), j=1..infinity)

Do you mean something like this?  You have a dataset representing points (x[i], y[i]) (approximately) on the graph of a rectangular function, say y = 1 for a < x < b, 0 otherwise, and you want to approximate (a+b)/2. 

Well, if the x[i] are in increasing order and y[i] = 1 for i1 <= i <= i2, 0 for i < i1 or i > i2, all you really know is that x[i1-1] <= a < x[i1] and x[i2] < b <= x[i2+1].  So in some sense the best available estimator for (a+b)/2 might be
(x[i1-1] + x[i1] + x[i2] + x[i2+1])/4.