For first question: 

f := x -> 36*x^6 + 2665*x^4 + 240*x - 675 + 4534*x^2 - 5836*x^3 - 516*x^5;

minimize(f(x), x = 0 .. 4, location);
#-675, {[{x = 0}, -675]}

evalf(maximize(f(x), x = 0 .. 4, location));
#703.9550742, {[{x = 3.800387934}, 703.9550742]}

Try:

ode := diff(U(z), z $ 4) + c^2*diff(U(z), z $ 2) + k*c*diff(U(z), z $ 2) - (3*U(z)^2 + a)*diff(U(z), z $ 2) = 0;
Order := 5;dsolve(ode, U(z), type = 'series');


#U(z) = U(0) + D(U)(0)*z + 1/2*(D@@2)(U)(0)*z^2 + 1/6*(D@@3)(U)(0)*z^3 + (U(0)^2*(D@@2)(U)(0)/8 - c^2*(D@@2)(U)(0)/24 - k*c*(D@@2)(U)(0)/24 + (D@@2)(U)(0)*a/24)*z^4 + O(z^5)

With initial conditions 

Order := 5;dsolve([ode, U(A) = A1, D(U)(A) = B1, (D@@2)(U)(A) = C1], U(z), type = 'series');

#U(z) = A1 + B1*(z - A) + 1/2*C1*(z - A)^2 + 1/6*(D@@3)(U)(A)*(z - A)^3 + (1/8*A1^2*C1 - 1/24*c^2*C1 - 1/24*k*c*C1 + 1/24*C1*a)*(z - A)^4 + O((z - A)^5)

As a workround using fourier transform:

(inttrans:-invfourier(int((inttrans:-fourier(sin(p*r), p, s) assuming (0 <= r))*sin(q*r)/(p*q), r = 0 .. infinity), s, p) assuming (q < p));

#-Pi*Dirac(p + q)/(2*p*q)

Try:

simplify(pdetest(sol, sys));

gives:

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Maybe this helps, see attached file:

General_formula_.mw

simplify(diff(int(JacobiSN(x, k)^2, x), x));

#JacobiSN(x, k)^2

See attached file:

EQ_v3.mw

Only solution,not  phase portrait.See atached file.

Solution.mw

One way is:

[seq(rhs(op(1, rootsq0[[n]])), n = 1 .. numelems([rootsq0]))];

Maybe like this:

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/integrals.mw .
 

Download integrals.mw

I don't have Maple 18.I don't understand what you mean by: h^k (0)=0 for k=0..n  ?

Solution by LaplaceTransform:

See attached file:
 

Download frac_eq.mw

frac_eq_2.mw

sol := dsolve((D@@2)(z)(t) + 2*D(z)(t) + z(t) = 2*exp(-t));
odetest(sol, (D@@2)(z)(t) + 2*D(z)(t) + z(t) = 2*exp(-t));
#0 ok.

In your case:

odetest(z(t) = t^2 + exp(-t), (D@@2)(z)(t) + 2*D(z)(t) + z(t) = 2*exp(-t));
#2 - 2*exp(-t) + 4*t + t^2

is not true.

For more info execute ?odetest.

genfunc:-rgf_pfrac(1/(z^2 + I)^2, z);

#(-1/4 + I/4)*sqrt(2)/(-sqrt(2)*I + sqrt(2) + 2*z) + I/(-sqrt(2)*I + sqrt(2) + 2*z)^2 + I/(sqrt(2)*I - sqrt(2) + 2*z)^2 + (1/4 - I/4)*sqrt(2)/(sqrt(2)*I - sqrt(2) + 2*z)

Using LaplaceTransform and Inverse LaplaceTransform. Fractional derivative is Riemann–Liouville sense.

restart;
v := t -> t;
plot([seq((inttrans:-invlaplace(s^alpha*inttrans:-laplace(v(t), t, s), s, t) assuming (0 < t)), alpha = -1 .. 1.5, 0.5)], t = 0 .. 10, view = [0 .. 9.5, 0 .. 5], legend = [seq('alpha' = alpha, alpha = -1 .. 1.5, 0.5)], color = [red, blue, green, yellow, cyan, magenta], axis[2] = [gridlines = [linestyle = dot]]);

Elzaki transform by Laplace Transform, see attached file.

Elziki_transforms_vs_2.mw