After posting I have noticed that Venkat already replied to Alex. 

@Axel Vogt As far as I understand, Venkat is using the upper limit as a parameter to define a ODE for my YFD(x) in terms of this new parameter he introduced.  The YDF limit is reached setting xi->infinity and setting the parameter x to whatever value we want to evaluate the YDF function at. The IC y(xi=x)=0 comes from the integral definition where the upper and lower limits are the same. Here is an worksheet that hopefully clarifes this. 

Download ParamYDF.mwParamYDF.mw

@ Venkat

Thanks so much for the optimization using the differential reprezentation to solve the integral. That really makes things faster! 

There is one problem I have noticed however. If you want to evaluate sol2 at high arguments. Certainly it is not well defined for x>100 as it stands right now, since in the parametric representation the upper limit of the integral was taken 100. However I have trouble estimating sol2 even at points much lower than this. 

@ Venkat, could you please also explain the role of:  if not type(evalf(x),'numeric') then
      'procname'(x) in your procedural definition of YFD. If that part is ommited we get back to the original error about xi. So what is the role of that line of code and how shold I interpret it? Thanks so much!

For your contributions!

@ Also, I dont really understand how to use initial x value as unknown parameter in the procedure output. Can u sketch a line of maple code of how to do that?

Indeed it works if I set Digits:=20 for stiff=true or rosenbrock so I can now use range and store the values. One question. I will extend this system to many more coupled ODEs. Say about 10000. However, I am interested in plotting only the first few, or combinations of the first few f[n,l]s, such as here f[0,0] and f[1,0]. In other words I wont need stored all the solutions for f[n,l] with arbitrary n and l in the range xini to xmax. With the option range the solution is stored, at least that is my understanding (all the f[n,l]s are computed and stored, so I might reach a problem with memory). So, how to I keep most of the dsolve solutions as procedures that are estimated once we need a their values for a certain x and only a few of them stored in memory for fast accessing when needed, such as the f[0,0](x), f[0,1](x) in the example above. 

@Preben Alsholm Yes, I am generating the ODEs using  two for nested loops, after n and after l. My functions are labeled f[n,l](x) and in each ODE the coefficients are mild functions of n or l (~linearly) but some of them have a very strong dependence on 1/x, such as 1/x^7. And i initialize at x=0.1, which I belive is the reason why I run into trouble with all the other methods. 

@Carl Love Thanks for the clarification.

@Carl Love Thanks, was not aware of this function. Yes, my question was actually related to the previous questions I have asked, since I have noticed Maple was keeping track of solutions that were as small as 10^-87, which for all intents and purposes were irrelevant. 

For all of your imput and help!

@acer 

You are the best!

Also, I have noticed that in your baseident you isolate the recurrsion relation you will be using. Is there a way to keep things general and just force Maple to use all identities available in FunctionAdviser when manipulating such expressions? Thanks!

Yes, I was thinking about your answer. However my level of Maple syntax is not that sophisticated yet, so I could not fully follow your example. Ill put here the expression that i need to simplify and perhaps you can help me with it?

(2l+2n-1)LaguerreL(n-2,l+1/2,y)-(2l+4n-2y-1)LaguerreL(n-1,l+1/2,y)

I can do this myself using the recurrsion relations, but I wanted Maple's help... How do i get Maple to use the FunctionAdviser in simplifying expressions?

And they are not assumed to obey l>0 and |m|≤l (as imposed by solving the Sturm-Liouville eigenvalue problem for the laplace operator on the unit sphere) when maple defines SphericalY. I believe that is the key to why my answer, already made good progress by working with the converted function:

However even if i can factor out now the l(l+1) eigenvalue, there is a trigonometric expresion of theta that does not simplify to 1. 

After some attempts the expression in question does simplify to 1 if one assumes we are not at the poles on the sphere, where laplace operator is ill defined. However would be nice if maple would just recognize the Laplace on Y_lm and spit out -l(l+1)Y_lm/r² without the need for manipulations.