making assumptions about x may eliminate some candidate solutions, you may also  remove rootofs with the AllSolutions option in solve, maybe something like

E[i, j] := Int(f[i](t)*u[j](t), t = -infinity .. infinity);
eqn[i, j] := u[j](x) = E[i, j];
k[i, j] := solve(eqn[i, j], x, AllSolutions) assuming x::positive;

also does the trick (but couldn't tell you which if any is to be commended)

subs(Result,alpha[2]);
                             -1.536774530

There must be a way to automate the following (with the use of the option discont=true), and someone else will certainly show how:

restart;
with(plots):
with(plottools):

fdiscont((x^2-4*x+5)/(x-2),x=1..3);

             [1.99972287449048180 .. 2.00074755004038485]

fdiscont helps you find the discontinuities of your function -- okay, so the asymptote is near 2 -- weird that Maple doesn't understand it's exactly 2.

Then I split the plot into two parts, so that Maple won't show a vertical line at the asymptote -- I then create my own asymptote by hand and give it the look I want (including dashes):

p1:=plot((x^2-4*x+5)/(x-2), x=1..1.99):
p2:=plot((x^2-4*x+5)/(x-2), x=2.001..3):
p3:=curve([[2,-10], [2,10]], color=red, linestyle=dash, thickness=2):
display({p1,p2,p3}, view=[1..3,-10..10]);

several methods discussed here:

http://www.mathworks.co.uk/matlabcentral/newsreader/view_thread/152947

assuming you meant to type * instead of = in the formula for H:

H:= (cos(n*L)*cosh(n*L))/(sin(n*L)-sinh(n*L)):
phi:=x->sin(n*x)-sinh(n*x)+H*(cos(n*x)-cosh(n*x)): 
int(phi(x), x=0..L);
  1/2 (4 exp(n L) sin(n L) - 4 exp(n L) sinh(n L)

         - 2 cos(n L) sin(n L) exp(n L)

         - cos(n L) cosh(n L) exp(2 n L)

         + 2 cos(n L) sinh(n L) exp(n L) + cos(n L) cosh(n L)

         - exp(2 n L) sin(n L) + exp(2 n L) sinh(n L) - sin(n L)

         + sinh(n L) + 2 cos(n L) cosh(n L) sin(n L) exp(n L))

        exp(-n L)/((sin(n L) - sinh(n L)) n)

I don't know much about this loglogplot, but if I run

loglogplot(10^x, x=1..2);

on my Windows/Maple13/Classic system I do get exponentials on the vertical axis.

You did not give us your code.

Explicit solutions are the exception rather than the rule. The following can be used to extract information from Maple's dsolve command.

restart;
with(DEtools);
infolevel[dsolve] := 4; 

dsys1 := {diff(u(t), t) = -(3/2-sqrt(2))*u(t)+5*v(t), diff(v(t), t) = -5*u(t)-(3/2+sqrt(2))*v(t)-20*cos(10*t)*w(t), diff(w(t), t) = -.5-3*w(t)+20*cos(10*t)*v(t)};

dsolve(dsys1);

-> Solving each unknown as a function of the next ones using the order: [u(t), w(t), v(t)]
-> Calling odsolve with the ODE, diff(diff(diff(y(x),x),x),x) = -1/4*(-48000*cos(10*x)^2*y(x)*sin(10*x)+7980*cos(10*x)^2+3100*diff(y(x),x)*sin(10*x)-728*diff(diff(y(x),x),x)*cos(10*x)+600*diff(diff(y(x),x),x)*sin(10*x)-319200*cos(10*x)^3*y(x)-19291*y(x)*cos(10*x)+4800*cos(10*x)^3*diff(y(x),x)+9090*y(x)*sin(10*x)+32000*cos(10*x)^2*diff(y(x),x)*sin(10*x)-960000*cos(10*x)*y(x)*sin(10*x)^2-1989*diff(y(x),x)*cos(10*x)+16000*sin(10*x)^2+(1600*cos(10*x)*sin(10*x)-96000*cos(10*x)^2*y(x)*sin(10*x)+80*diff(diff(y(x),x),x)*sin(10*x)+240*diff(y(x),x)*sin(10*x)+2020*y(x)*sin(10*x)+3200*cos(10*x)^3*diff(y(x),x)+48*diff(diff(y(x),x),x)*cos(10*x)+274*diff(y(x),x)*cos(10*x)+606*y(x)*cos(10*x))*2^(1/2))/(3*cos(10*x)+20*sin(10*x)+2*2^(1/2)*cos(10*x)), y(x), singsol = none
Methods for third order ODEs:
--- Trying classification methods ---
trying a quadrature
trying high order exact linear fully integrable
trying differential order: 3; linear nonhomogeneous with symmetry [0,1]
trying high order linear exact nonhomogeneous
trying differential order: 3; missing the dependent variable
-> pFq: Equivalence to the 3F2 or one of its 3 confluent cases under 'a power @ Moebius'
trying a solution in terms of MeijerG functions
-> Try computing a Rational Normal Form for the given ODE...
<- unable to resolve the Equivalence to a Rational Normal Form
checking substitution methods
trying differential order: 3; missing the dependent variable
-> pFq: Equivalence to the 3F2 or one of its 3 confluent cases under 'a power @ Moebius'
trying a solution in terms of MeijerG functions
-> Try computing a Rational Normal Form for the given ODE...
<- unable to resolve the Equivalence to a Rational Normal Form
checking substitution methods
Removed one radical, (-4545-979291*t^5-979291*t+91455*t^2-1010*2^(1/2)+190990*2^(1/2)*t^2+606*2^(1/2)*t+1212*2^(1/2)*t^3+606*2^(1/2)*t^5-91455*t^4+4545*t^6+604618*t^3-190990*2^(1/2)*t^4+1010*2^(1/2)*t^6)*u+_u[1]*(-21945*t+54165*t^3+102165*t^5+24000*t^9-421900*t^2-169200*t^4+246300*t^6+41350*t^8+30000*t^10+25370*2^(1/2)*t^7+8400*2^(1/2)*t^6+92110*2^(1/2)*t^5+7200*2^(1/2)*t^4+76110*2^(1/2)*t^3+1200*2^(1/2)*t^2+1370*2^(1/2)*t+3000*2^(1/2)*t^8-17750+50055*t^7+8000*t^9*2^(1/2)-600*2^(1/2))+_u[2]*(-48200*t-13500*t^2+21000*t^4+69000*t^6+58500*t^8+27000*2^(1/2)*t^8+4800*2^(1/2)*t^7+38000*2^(1/2)*t^6+7200*2^(1/2)*t^5+4800*2^(1/2)*t^3+3000*2^(1/2)*t^2+1200*2^(1/2)*t+22000*2^(1/2)*t^4-7500-162800*t^3-169200*t^5-12800*t^7+71800*t^9+30000*t^11+16500*t^10+1200*t^9*2^(1/2)+7000*t^10*2^(1/2)-1000*2^(1/2))+_u[3]*(-20000*t^2-25000*t^4+25000*t^8+20000*t^10+5000*t^12+5000*t^9*2^(1/2)+10000*2^(1/2)*t^5+1000*2^(1/2)*t-5000+10000*2^(1/2)*t^7+5000*2^(1/2)*t^3+1000*t^11*2^(1/2)+7500*t^3+15000*t^7+1500*t^11+1500*t+15000*t^5+7500*t^9)
trying differential order: 3; missing the dependent variable
checking if the LODE is of Euler type
trying Louvillian solutions for 3rd order ODEs, imprimitive case
-> pFq: Equivalence to the 3F2 or one of its 3 confluent cases under 'a power @ Moebius'

Maple then gives you an equivalent representation of your system in terms of one third-degree ODE in v(t). Run the above code and you can see that it's unwieldy.

I don't get the memory allocation error problem.

I am guessing that EU must be something like: Px*U(x)+Py*U(y), where U() is your utility function.

if your utility function is quadratic, since your constraint linear, you should/must be able to solve the problem in closed form.

if U(x)=1-x^2, say...

Takes a short while to produce but looks pretty smooth:

contourplot( max(min(-z*(1/(x^2+z^2)^(3/2)-1), 1000),-1000), 
x=-0.3..0.3, z=-0.3..0.3, numpoints=50^3,axes=boxed);

The max min trick is pretty cool. Its usefulness is particularly apparent, in this example, for small values of numpoints:

plot3d( max(min(-z*(1/(x^2+z^2)^(3/2)-1), 1000),-1000),x=-0.3..0.3,z=-0.3..0.3, 
numpoints=50,axes=boxed,lightmodel=light2,orientation=[27, 54]);

versus

plot3d( -z*(1/(x^2+z^2)^(3/2)-1), x=-0.3..0.3, z=-0.3..0.3,
view=[DEFAULT,DEFAULT,-1000..1000], numpoints=50,
axes=boxed,lightmodel=light2,orientation=[27, 54]);

Okay, I confess, I'm really practicing my image-posting skill.

you're absolutely right, I copied the wrong line of code, apologies.

I think. You can see that by gradually raising the range of t from, say, 0..1 to 0..4, the following sums it up:

plots[odeplot](sol,[t,Theta1(t)],0..2);
plots[odeplot](sol,[t,Theta2(t)],0..2);

in this case, solve does it:

y=(-4+4*d)*x^2+(4-4*d)*x-3*d+2;
solve(%, x);

one way to do it: numpoints

the default numpoints is 50, increase it for better effect:

contourplot(-z*(1/(x^2+z^2)^(3/2)-1),x=-3..3,z=-3..3,grid=[50,50],contours=10,numpoints=100);

with(plots):
f := x -> x^(-2):
logplot(f(x), x=0..2);