In Maple 11 the result of dsolve(ode) is

y(x) = _C1*LegendreP(j,k,x)+_C2*LegendreQ(j,k,x)

Your initial conditions are at x=1, but there's a problem: the differential

equation is singular at x=1, with a regular singular point there. 

 

> DEtools[indicialeq](ode, x,1,y(x));

x^2-1/4*k^2 = 0

 

(Hmmm, looks like a problem in MaplePrimes here: this should be x^2-1/4*k^2=0) 

So the indicial roots are (+/-) k/2.  Assuming k > 0, most solutions will have |y| -> infinity as x -> 1, and those that don't (which I believe are the ones with _C2 = 0) will have |y| -> 0 as x -> 1.  There are no solutions with y -> 1 as x -> 1.

Actually this is a case of integration by parts.

> with(student):

  J := int(int(f(s),s=0..t),t=0..x);

 intparts(J, op(1,J));

int(f(s),s = 0 .. x)*x-Int(f(t)*t,t = 0 .. x)

> factor(combine(subs(s=t,value(%))));

int(f(t)*(-t+x),t = 0 .. x)

 

 

You could try this:

> S:= expand([seq](``(j-1<x and x <= j, (j-x)*`if`(j=1,0,M[i,j-1])+(x-j+1)*M[i,j]), j=1..18));

And then to get your 10 functions

> f:= table([seq](unapply(piecewise(op(S)), x), i=1..10));

 

 

I don't know why my suggestion wouldn't work for you.  Also, if you got a complaint about the final value in the for loop, what you entered couldn't have been

for i from 1 to 100 do

Perhaps if you upload your worksheet we could see what's wrong.

 

To go to a new line without triggering execution, press Shift and Enter at the same time. 

If you were trying to do multi-line loops in 1D Maple input without Shift+Enter, you must have seen

Warning, premature end of input, use <Shift> + <Enter> to avoid this message.

 

 

> type(L, list(positive));

Yes, that should work fine (assuming evalf is able to give a numerical value to each entry of M).

Yes, it's a bug.  The work-around is to do it in Classic.

First, let's find the eigenvalues, and do what simplification can be done with the benefit of the side relation a^2 + b^2 = 1.

> Evals := simplify(Eigenvalues(U), {a^2 + b^2 = 1});

Evals := RTABLE(166792836,MATRIX([[1/2*I+1/4*(-4+(8-8*a)^(1/2))^(1/2)], [1/2*I-1/4*(-4+(8-8*a)^(1/2))^(1/2)], [1/2*I+1/4*(-4-(8-8*a)^(1/2))^(1/2)], [1/2*I-1/4*(-4-(8-8*a)^(1/2))^(1/2)]]),Vector[column]) So far so good: you might note that (except for some special values of a) these are all distinct, and they are roots of the characteristic polynomial of U whenever a^2 + b^2 = 1:

> P := CharacteristicPolynomial(U,lambda):
> map(t -> simplify(eval(P, lambda = t),{a^2+b^2=1}), Evals);

RTABLE(149986320,MATRIX([[0], [0], [0], [0]]),Vector[column]) Now, since Eigenvectors(U) doesn't produce a very nice result, here's a trick I sometimes use: for each eigenvalue lambda, use LinearSolve on the first 3 rows of U - lambda*Identity.

> U1:= U - Evals[1]*IdentityMatrix(4):
  V1:= LinearSolve(U1[1..3,1..4],<0,0,0>);

V1 := RTABLE(150233796,MATRIX([[_t[1]], [1/(a*(-4+2*(-2*a+2)^(1/2))^(1/2)-(-4+2*(-2*a+2)^(1/2))^(1/2) +2*I*b-b*(-4+2*(-2*a+2)^(1/2))^(1/2))*(2*I*a+a*(-4+2* (-2*a+2)^(1/2))^(1/2)-2*I-(-4+2*(-2*a+2)^(1/2))^(1/2) +b*(-4+2*(-2*a+2)^(1/2))^(1/2)+2*I*(-2*a+2)^(1/2))*_t[1]], [1/4*I*(2*I+(-4+2*(-2*a+2)^(1/2))^(1/2))*_t[1]], [1/4*I*(2*I-(-4+2*(-2*a+2)^(1/2))^(1/2))/(a*(-4+2* (-2*a+2)^(1/2))^(1/2)-(-4+2*(-2*a+2)^(1/2))^(1/2)+ 2*I*b-b*(-4+2*(-2*a+2)^(1/2))^(1/2))*(2*I*a+a*(-4+2* (-2*a+2)^(1/2))^(1/2)-2*I-(-4+2*(-2*a+2)^(1/2))^(1/2) +b*(-4+2*(-2*a+2)^(1/2))^(1/2)+2*I*(-2*a+2)^(1/2))*_t[1]]]), Vector[column]) The solution contains an arbitrary parameter _t[1]: set it to 1.

> V1:= eval(V1,_t[1]=1);

Check that this is an eigenvector for Evals[1]:

> map(simplify,U.V1-Evals[1]*V1,{a^2+b^2=1});
  map(normal,%);

RTABLE(165361656,MATRIX([[0], [0], [0], [0]]),Vector[column]) You can get the other eigenvectors similarly.

As far as I know, the gridlines option is only for 2D plots. However, you can make your own grid lines "by hand":

> with(plots):
  P1:= plot3d(x^2-y^2, x = -1 .. 1, y = -1 .. 1):
  Grid:= seq(seq(spacecurve([[i/2,j/2,-1],[i/2,j/2,1]], 
     colour=gray),i=-2..2),j=-2..2),
         seq(seq(spacecurve([[i/2,-1,j/2],[i/2,1,j/2]], 
     colour=gray),i=-2..2),j=-2..2),
         seq(seq(spacecurve([[-1,i/2,j/2],[1,i/2,j/2]], 
     colour=gray),i=-2..2),j=-2..2):
  display([Grid,P1],axes=boxed);

For a 3D bar chart, you might try matrixplot in the plots package with the option heights=histogram.

Hmmm: my copy of Goldstein's "Classical Mechanics" doesn't have this. I guess that's because it was printed in 1965... Well, I think the place to start off would be the tensor package. In particular, you might look at the examples for the "Einstein" command in that package.

> solve({ your equations });

Since the characteristic polynomial is a 5th degree polynomial in H, you can't just "get rid of" the RootOf: it's the best way to represent the solution. Besides what John did, you could try these approaches:

> P := CharacteristicPolynomial(H, z):
  implicitplot(P, F = 0 .. 0.1, z = -1 .. 1,   
     gridrefine=3, crossingrefine=3);
  plot([seq](RootOf(P,z,index=j),j=1..5),F=0..0.1,
     discont=true);

Note that except for F very close to 0, the curves seem to be very well approximated by straight lines. Those may be obtained as follows:

> A := asympt(eval(P,z=p+q*F), F);

A := (q^5+1858126385941085927943/2980232238769531250*q -1012333438017/12207031250*q^3)*F^5+ (254288177326787325327/5960464477539062500- 1154456949273/48828125000*q^2+5*q^4*p+5/12*q^4 -3037000314051/12207031250*q^2*p+ 1858126385941085927943/2980232238769531250*p)*F^4+ (5/3*q^3*p-1154456949273/24414062500*q*p+10*q^3*p^2 -5120339292773/2343750000000*q-3037000314051/ 12207031250*q*p^2+115/1728*q^3)*F^3+(-1012333438017/ 12207031250*p^3+115/576*q^2*p+10*q^2*p^3+5/2*q^2*p^2 -5120339292773/2343750000000*p-1154456949273/ 48828125000*p^2+475/93312*q^2-7925066173943/ 126562500000000)*F^2+(5/3*q*p^3+5*q*p^4+35/186624*q+ 475/46656*q*p+115/576*q*p^2)*F+5/12*p^4+475/93312*p^2 +p^5+35/186624*p+115/1728*p^3+1/373248

> qs:= [solve(coeff(A,F,5),q)]:
  ps:= [seq](solve(coeff(A,F,4),p), q = qs):
  Es:= evalf(zip((q,p) -> p+q*F,qs,ps));

Es := [-.6842596372e-1, -8.635346757*F-.6041065906e-1, 8.635346757*F-.6041065906e-1, -.1137096922-2.891563868*F, -.1137096922+2.891563868*F]

Assuming you are using the old-style "matrix" rather than the new-style "Matrix":

> A := copy(B);

If you were using "Matrix":

> with(LinearAlgebra):
  A := Copy(B);

For plotting, try DEplot in the DEtools package. See the first example in the help page ?DEplot. For finding the equilibrium solutions, I suggest you review what an equilibrium solution is. You don't need Maple for that in this example.