@Axel Vogt Thanks! I will try some more tricks with combine(), Re() etc,. Also, the problem with elminating k1, k2 is that this is not possible for all 603 integrals. I've been playing around with exporting maple to MathML, since mathematica can read in those files. Unfortunetly, I can't find a simple command to export a MathML expession to a .mml file. I tried writeto() which saves as a .txt file? And then tried to read in with MMa's Import[%, "mathML"] with no luck.  

Edit: I can convert to Mathematica pretty easily now. I wrote a simple script to replace all "sin(x)" with "Sin[x]", etc....  

@Markiyan Hirnyk yes, thats what I got too. Perhaps I made a typo when converting from Mma to Maple. 

@ecterrab awesome! thanks for the quick reply.

@ecterrab Thanks so much! You saved me a world of pain in trying to figure this out!

seq is much much faster than for loops. You should try to re-write the code using only seq. If you have a multi-core computer, you could try Threads:-Seq or Grid:-Seq. 

@ecterrab Great thanks for the quick reply! 

@Markiyan Hirnyk So in conclusion there is no way to save memory on this since interpolate has been removed for IVP's?

btw this was related to an earlier post http://www.mapleprimes.com/questions/149602-Solving-First-Order-Linear-Coupled-DE-Matrix-Form where I was trying to solve  a set of first order, linear ODE's in matrix form, I would love to be able to solve all 961 equations for more accuracy, and the minimun is aroudn 132, I can od about 170 before I get a memory error. 

X:=Vector([seq(x[i](t),i=1..(L+1)^2)]);
diff~(X,t)=mat.X; #Your system
eqs:=convert(Equate(lhs(%),rhs(%)),set);
t0:=1: #Initial t
inits:={seq(x[i](t0)=x0[i],i=1..(L+1)^2)};

res:=dsolve(eqs union inits,numeric,parameters=[seq(x0[i],i=1..(L+1)^2)],range=1..15);
res(parameters=[1,2,3,.1,.2,.3,4,5,6]);
plots:-odeplot(res,[Re(x[1](t)),Im(x[1](t))],1..15,refine=1

From http://www.maplesoft.com/support/help/Maple/view.aspx?path=dsolve/rkf45

Also the matrix DE is a system of 961 (testing on 161 currently)  linear, coupled, complex valued DEs

that works surprisingly fast! Thanks! But with the parameters option, can I pass complex numbers? i.e.

coef := x0[1] = .1126537-0.2275632e-2*I, x0[2] = 0.7460254e-1+0.6902865e-1*I, x0[3] = .1339367-1.767657*10^(-7)*I, x0[4] = -0.6367112e-1+0.6881254e-1*I, x0[5] = 0.4677631e-2+0.8336431e-1*I, x0[6] = .1101614+.1061479*I, x0[7] = 0.5578725e-1+0.5915829e-2*I, x0[8] = -0.9790946e-1+.1016116*I, x0[9] = -0.3979998e-2-0.7092869e-1*I

c:=convert(%,list);

res:=dsolve(eqs union inits,numeric,parameters=[seq(x0[i],i=1..(L+1)^2)],range=1..15);
res(parameters=[c] );

that works surprisingly fast! Thanks! But with the parameters option, can I pass complex numbers? i.e.

coef := x0[1] = .1126537-0.2275632e-2*I, x0[2] = 0.7460254e-1+0.6902865e-1*I, x0[3] = .1339367-1.767657*10^(-7)*I, x0[4] = -0.6367112e-1+0.6881254e-1*I, x0[5] = 0.4677631e-2+0.8336431e-1*I, x0[6] = .1101614+.1061479*I, x0[7] = 0.5578725e-1+0.5915829e-2*I, x0[8] = -0.9790946e-1+.1016116*I, x0[9] = -0.3979998e-2-0.7092869e-1*I

c:=convert(%,list);

res:=dsolve(eqs union inits,numeric,parameters=[seq(x0[i],i=1..(L+1)^2)],range=1..15);
res(parameters=[c] );

@mehdi_mech  Yes, if we treat y'=y'(x) its circular to find x in terms of y'(x)  [ y'(x(y'(x....))) doesnt make sense]. But if you merely want everything in terms of a new independant variable, y', which happens to be dy/dx in the previous variable system then we can assign x(y') and y(y') that satisfy the above post it amounts to just a change of variables. I don't think this serves and practical reason but we can obtain x(y') and y(y') up to a numerical constant from y' alone. 

@mehdi_mech  Yes, if we treat y'=y'(x) its circular to find x in terms of y'(x)  [ y'(x(y'(x....))) doesnt make sense]. But if you merely want everything in terms of a new independant variable, y', which happens to be dy/dx in the previous variable system then we can assign x(y') and y(y') that satisfy the above post it amounts to just a change of variables. I don't think this serves and practical reason but we can obtain x(y') and y(y') up to a numerical constant from y' alone. 

@Markiyan Hirnyk  given that we found y(y') and y''(y') we can find x(y') from  int( 1/y'' dy' ), but I think my code above also treats x as x(y') since

df/dy' = df/dx * dx/dy  = df/dx *1/y', and we can rearrange and integrate for x(y').

We could do more work and find the inverse function without calculus, if we were given a specific function. i.e. y'=sin(x), x=arcsin(y') is the inverse. So I tentatively say yes, it's  possible and already implemented. 

@Markiyan Hirnyk  given that we found y(y') and y''(y') we can find x(y') from  int( 1/y'' dy' ), but I think my code above also treats x as x(y') since

df/dy' = df/dx * dx/dy  = df/dx *1/y', and we can rearrange and integrate for x(y').

We could do more work and find the inverse function without calculus, if we were given a specific function. i.e. y'=sin(x), x=arcsin(y') is the inverse. So I tentatively say yes, it's  possible and already implemented. 

@Axel Vogt I was pointing out another mistake; sometimes Maple is off by a negitive sign like the (3,-3),(3,3) case and other times it is off mixing up real/img. parts like (0,0),(3,3). 

I proved that the imaginary part must be zero...not the real part, can you point out where I might have made a mistake? 

What other conventions are there? I'm using the standard physics one found in Griffiths, Mathematica, Wikipedia, etc. 

I wrote my own SphericalY proc that fixes the issue (I think) and heres the file that compares their results in the (3,-3),(0,0) case. I noticed they only disagree when I have something of the form Y(L,-L,theta,phi).

bug_in_sphericals.mw