Note that the right side of your equation is proportional to x, so your solution must be
y(x) = c*x for some c. So just plug that in:

> eq:= y(x) = 7/8 * x + 1/2 * int (x*t*y(t)^2,t=0..1);
   eval(eq, y = (x -> c*x));

And then solve for c:

> solve(%, c);

1, 7

Probably the equation you're really interested in is somewhat more complicated, though.
In your proposed method of solution:

> H:= (u,p) -> (1-p)*(7/8 * x - u(x)) 
       + p*(7/8 * x + 1/2*int(x*t*u(t)^2, t=0..1) - u(x));
   U:= unapply(add(u[i](x)*p^i, i=0..5), x);
   S:=series(H(U,p),p,6);
   for i from 0 to 5 do u[i]:= unapply(solve(coeff(S,p,i),u[i](x)),x) end do;

And then the (approximate) solution to your equation:

> eval(U(x), p=1);

with(plots):

animate(pointplot3d,[[seq(seq(seq([r*t,i*Pi,j*Pi],r=0..1,0.1),i=0..2,1/5),j=0..1,1/10)], coords=spherical], t = 0 .. 10,symbol=diamond,shading=zhue);

You might look at the paper of Boswell and Glasser "Solvable Sextic Equations", http://www.arxiv.com/abs/math-ph/0504001

They mention using Maple; I wonder if they have Maple code available.

plot3d is for surfaces, not curves.  spacecurve in the plots package is for plotting curves in 3D.  You could try something like this:

> P1 := 1 + 2*Q1^2; P11:= 2 + Q1^2;
   P2:= 3 + 4*Q2^2; P22:= 5 + 3*Q2^2;
   with(plots):
   solve(P1 = P11);
   m11, m12:= (min,max)(%);
   solve(P2 = P22);
   m21, m22:= (min,max)(%);
   display([spacecurve([Q1, 0, P1], Q1 = -3 .. m11, colour=red, linestyle=dash), 
                 spacecurve([Q1, 0, P1], Q1 = m11 .. m12, colour=red),   
                 spacecurve([Q1, 0, P1], Q1 = m12 .. 3, colour=red, linestyle=dash),
                 spacecurve([0, Q2, P2], Q2 = -3 .. m21, colour=blue),
                 spacecurve([0, Q2, P2], Q2 = m21 .. m22, colour=blue, linestyle=dash),
                 spacecurve([0, Q2, P2], Q2 = m22 .. 3, colour=blue)]);

You mean like this?

> plots[pointplot3d]([seq(seq(seq([r,i*Pi,j*Pi], r = 0 .. 1, 0.1), i=0..2, 1/5), j = 0 .. 1, 1/10)], coords=spherical, symbol=diamond, shading=zhue);

For your solvable quintic, solve works:

solve(32*c^5-16*c^4-32*c^3+12*c^2+6*c-1, explicit);

For the sextic, it doesn't. 

Square brackets [] are for subscripts and lists, never for grouping.  Braces {} are for sets, not for grouping.  Defining a[j] defines only (literally) a[j], not a[anything else].

I think you want something like this:

> r:= 0.05; Deltat:= T/N; DeltaS:=S/M; T:=0.4; q:=0; K:=1.1; S:=2; M:= 10; N:=10; 
   sigma:= 0.1;
   for j from 0 to M do f[N,j]:= max(K - j*DeltaS,0) end do;
   for j from 0 to M do 
     a[j]:= (-1/2 *(r-q)*j*Deltat + 1/2*sigma^2*j^2*Deltat)/1 + r*Deltat;
     b[j]:= (1 - sigma^2 * j^2 * Deltat) / 1 + r * Deltat;
     c[j]:= (1/2 * (r-q) * j * Deltat + 1/2 * sigma^2 * j^2 * Deltat) / 1 + r * Deltat
   end do;
   for i from N-1 to 0 by -1 do
     for j from 1 to M-1 do
       f[i,j]:= a[j]*f[i+1,j-1] + b[j]*f[i+1,j] + c[j]*f[i+1,j+1]
     end do;
     f[i,0]:= b[0]*f[i+1,0] + c[0]*f[i+1,1]; 
     f[i,M]:= a[M]*f[i+1,M-1] + b[M]*f[i+1,M];
 end do:

Something like this?

> with(plots):
P1:= pointplot([[2,3],[-3,1],[-1.5,-2.5],[0,0]],colour=[green,red,blue,magenta],
    symbol=solidcircle,symbolsize=8):
P2:= plot([[[0,3],[2,3],[2,0]],[[-3,0],[-3,1],[0,1]],[[-1.5,0],[-1.5,-2.5],[0,-2.5]]],
    linestyle=dash,colour=[green,red,blue]):
P3:= textplot([[2.2,3.3,"(2,3)"],[-2.8,1.3,"(-3,1)"],[-1.8,-2.8,"(-1.5,-2.5)"],   
    [0.4,0.4,"(0,0)"]],colour=[green,red,blue,magenta]):
display(P1,P2,P3,scaling=constrained,view=[-3.9..3.9,-3.9..3.9],labels=[x,y]);

You could just as easily use plot.

> plot([seq([seq([0.1e-1*2^(k-1),F[i,k]], k=dmin .. dmax)], i=1..3)], 
      colour=black, linestyle=[solid, dash, dot]);

Knowing that sqrt(5) is involved in the factoring:

> X1:= (s^2+2)*s/(1+s^4+3*s^2);
   factor(X1,{sqrt(5)});
   convert(%, parfrac);
   inttrans[invlaplace](%,s,t);

                         1/2                             1/2
                      t 5            1/2              t 5
         cos(t/2) cos(------) + 1/5 5    sin(t/2) sin(------)
                        2                               2

Alternatively, you could factor the denominator over the (floating-point) reals.

> numer(X1)/factor(denom(X1),real);
   inttrans[invlaplace](%,s,t);

.7236067977*cos(.6180339888*t)+.2763932022*cos(1.618033989*t)

> identify(%);

 /       1/2\     / 1/2      \      /       1/2\     /       1/2\
 |      5   |     |5         |      |      5   |     |      5   |
 |1/2 + ----| cos(|---- - 1/2| t) + |1/2 - ----| cos(|1/2 + ----| t)
 \       10 /     \ 2        /      \       10 /     \       2  /

There's not much that will stop the classical rk4 method, except the differential equation becoming undefined or the solution becoming complex.  As Preben says, we'd need the actual system to tell more.

What you have is not a polynomial.  I think what you want is this:

> f := mul(1/(1-x^j), j=1..10);
   series(f, x, 11);

By definition, permutations can contain the same element only once.  What you are asking for are elements of P^3 (i.e. all ordered triples of elements of P). They can be obtained using cartprod in the combinat package, or simply using a nested seq.

> with(combinat):
   Iterator:= cartprod([P,P,P]):
   [seq(Iterator[nextvalue](), i=1..nops(P)^3)];

or

  [seq(seq(seq([a,b,c],a=P),b=P),c=P)];

You're probably assuming sigma > 0 (or at least sigma real and nonzero), but Maple makes no such assumption unless you tell it to do so.  Without such an assumption, the integral might not converge.

> int(1/8*2^(1/2)*exp(-1/2*x^2/sigma^2)/(Pi^(5/2)*sigma)+1/8*2^(1/2)*exp(-1/2*x^2/sigma^2)*Q*exp(I*k*x)*exp(I*k*theta)/(Pi^(3/2)*sigma)+1/8*2^(1/2)*exp(-1/2*x^2/sigma^2)*Q/(Pi^(3/2)*sigma*exp(I*k*x)*exp(I*k*theta)), x = (-infinity .. infinity)) assuming sigma > 0;

Something like this, perhaps, for the matrix

[1/2  1 ]
[1/2  1 ]

>  with(plots):
     P1:= plot({[[-1.2,-2.4],[1.2,2.4]],[[-2.4,1.2],[2.4,-1.2]]}, colour=black):
     P2:= plot({[[3.6,2.4],[8.4,-2.4]],[[3.6,-2.4],[8.4,2.4]]}, colour=black):
     P3:= plot([[-1,1/2],[-1/4,2],[3/4,3/2]], colour=black, linestyle=dash):
     P4:= plot({[[-1,1/2],[6,0]],[[-1/4,2],[6+15/8,15/8],[3/4,3/2]]}, colour=red, 
thickness=2):
     L:= textplot([[-1.8,.5,N(A)],[-.1,-1.5,C(A^T)],[5,-1.7,C(A)],[7.5,-2.2,N(A^T)]]):
     display(P1,P2,P3,P4,L, axes=none, scaling=constrained);