@sursumCorda OP said multivariable polynomial, and I think it is OK for that. I suppose it would be more robust if you expanded p first, but the meaning of the answer depends on the order of the terms that is given, so I assumed that it was given as a sum of terms.

@SHIVAS Ans[k] is the output from pdsolve, which is a module - it can't be used in the same way as the output from dsolve. See the help page ?pdsolve,numeric for how to use the module to create plots and find values.

@Tamour_Zubair 

solve(subs(y(t)=y,sol),t); gives t as a function of y.

@Tamour_Zubair There isn't an analytical solution for y as a function of t; the RootOf is the best you can do. On the other hand, from the implicit solution you can get t as an analytical function of y. 

Please check your boundary condition D(f)(0) = lambda + xi*(D^2)(f)(0). Note that if you wanted the second derivative and not the square of the first derivative, then you need D(f)(0) = lambda + xi*(D@@2)(f)(0)

@NeraSnow Your procedure to upload the worksheet was correct, but the website is not working correctly today, as happens sometimes.

@Carl Love Yes, I tried something like this expecting the new inline loops might be more efficient, but was disappointed (as I've found a couple of times before).

There is also the issue of counts>9, which I missed originally. So "1111111111111" gives "131" in my code, but "11111111111111" in yours.

@sursumCorda Arrays are usually efficient, but here adding on to the end with the ,= operation is part of the secret.

@AngieS7 It means it can't find a solution. You can help it by providing ranges (or an initial guess), but it might mean there is no solution

@AngieS7 You can tell fsolve to only look for solutions in specific ranges by using something like

fsolve({eq1, eq2, eq3, eq4, eq5, eq6}, {A=0..5, B=3..7, C=0..infinity, omega=0..10, x1=0..1, x2=0..1});

@AngieS7 It is hard to tell why this hasn't worked without you uploading your worksheet. You can do that with the green up-arrow in the Mapleprimes editor.

Edit: looks like you might have changed vb:=0.5 to vb:=5

@hojat Lambda still has a 1[1] and a 2[1].

@hojat Please confirm that you are minimizing a function of 12 variables (from indets(Lambda,name) ). If so this will be an extremely difficult problem and I'm guessing it generates some complex numbers because it doesn't know enough about these values. What are the approximate values of these variables? For nonlinear problems, an initialpoint is almost certainly needed.

@ianmccr This recasting into alternate forms happens in several situations, including just printing the code of a procedure with a semicolon. Since these are often less readable, I find this frustrating but don't know of a workaround.

@sursumCorda I see there that the symmetric NAG routine is recommended over the general one because it is "more accurate", so a NAG general routine is probably not useful here (can't increase digits then) though the specialized ones may be. Maple probably has a special routine for exp (I didn't check), which is why @Axel Vogt's method works.

Detecting the non-diagonalizable case will always be difficult if Jordan forms are used, because that is a difficult computation. Eigenvectors returns a "zero eigenvector" in that case but how close to zero does it have to be to be detected? [Edit: I see Eigenvectors can just output the true eigenvectors, so has some detection method for this.] These situations, including when two numerical eigenvalues are the same, rely on some careful assessments and are tricky in general - maybe the linalg method that worked well in this case would not work so well in other cases.

However, the fact that it is so bad in quite a few cases means it needs some improvement for sure.