@NIMA112 That's why it is a reply and not an answer. As I said "I don't know precisely what you want, since you do not show any calculation that didn't work." Your title mentions a numerical error but your worksheet does not show any explicit error. What is your question? What exactly do you think the "singularity problem" is? I am guessing that the physical problem might have a singularity and therefore no amount of mathematics will make it go away.

I don't know precisely what you want, since you do not show any calculation that didn't work. But visually there does seem to be a singularity in the physical problem at (-1,1.5) where the contour lines converge. So I guess you can go closer by using a higher setting of Digits. Note that you have mispelled Digits so it is not actually at 30.

@Nicole Sharp I don't use the Maple add-on in Excel so don't have an answer for all the inconsistencies you see. So just two comments that might be helpful

Since Excel works in hardware precision (about 15 digits), if Maple passes it a result with 32 digit precision Excel must truncate that in further calculations, so how the extra ones are displayed or if they are just removed is a moot point.

9*10^9 is exact and not a floating point number with finite precision, so displays as "9000000000". evalf converts it to a floating point number. But 9.*10^9 has an explicit decimal point and so is a floating point number. In Maple, numbers created with the "e" or "E" exponent notation are always floating point, even if there is no explicit decimal point.

@SHIVAS As in my last reply, pdsolve can't find an exact solution, so I don't know how to find Theta_exact. Do you have some reason to believe there is an analytical solution? If so, you can perhaps use hints or other methods to help pdsolve it, but you would need to know or guess something about the solution.

@SHIVAS 

For your new equation,

pdsolve(OdeSys, Theta(X, tau))

gives no solution so it is unlikely that one exists.

@2cUniverse 

I put almost all the code in the startup code region, and in the click, changed etc code regions I just call a procedure defined in the startup region. This means you aren't debugging multiple bits of code in multple locations, and you can have extensive code without any problem.

I loaded the packages in the startup region and it seemed to work consistently.

pb-App_0.1.mw

@mmcdara Interesting approach. In the small amount of reading I did, I recall something about the fact that values need not be near the corresponding eigenvalue. For contour plots, I'm not sure the results are invalid, just that you might need very fine countours to catch the unusual locations?

@rzarouf These sorts of full evaluation problems for making plots are common and are best solved by using the procedure form of the plot call, though you also found a bug in Norm, for which I will submit an SCR.

Download pseudospectrum2.mw

@Carl Love Thanks. Sometimes I hope the OPs provide at least some information to help. (If it's not in Horn and Johnson it seems it must be a bit exotic.) Turns out it was easier than I thought to calculate.

Please explain what a pseudospectrum is. If you just mean to plot the eigenvalues of a matrix in the complex plane, that can easily be done with complexplot.

@SHIVAS You need to solve and create the plots in a loop and then combine them with plots:-display. For example if your plots for different Nr are eta[1],eta[2] etc, then they are combined by something like plots:-display([seq(eta[i],i=1..3)])

@treverona Once you have a symbolic power like your d then Maple doesn't consider it as a polynomial, so the usual commands to find degrees don't work.

@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.