Here it is too_many_levels_of_recursion.mw

Thanks @Scot Gould . I've been looking for a resource whoch explains how to implement QM in MAPLE. 

I'm currently in the process of playing around with anharmonic potentials and obtaining eigen energies and wavefunctions. Your intial conditions work much better than my own.

My problem is a bit tricky as it involves Heun functions which are implemented in MAPLE. 

Any assistance such as the one you've just given is hugely appreciated :-)

@Preben Alsholm Thank you. Really appreciate it.

@Kitonum Sorry to disturb but was just wondering if you've given this problem a go. It's fine if you haven't and I really appreciate what you have already done

@AmusingYeti Sorry I am a real novice at this and really have no idea how it's done. Could you please add your brilliant suggestion to the code and give us the final code. Thanks in advance

@Preben Alsholm Could you please give me the full code(by adding your bit to the code I already have thanks to Kitonium). Really sorry for causing trouble but I'm really a beginner and these things often seem impossible to me. I simply can't seem to get them normalized.Thanks in advance:-)

@Kitonum I'm new here and didn't realise that iterrupting a thread isn't recommended. Also sorry for plagiarising you. 

Normalizing means that the square of the wavefunction integrates to 1 between -infinty and infinity(see amusing yeti post below).

Thanks and sorry again

@Kitonum Thank you very much. I am very much a beginner and really appreciate any help I can get.

Do you mind if I ask another quick one. Say just E=1.06538 has a different initial conditon(i.e D(psi)(3)=-1 instead) is it still possible to plot on same graph. If so how do I change the code?

Thanks in advance.

@Kitonum Ya I know it does mate. 

The problem is rather more trivial. Suppose I change the energy to 1.06538,5.74795,8.35258 I get different wavefunctions. I want all these wavefunctions on the same plot so I can compare their shapes 

Also I need them to be normalized. Thanks

Sorry the above code looks messy. Here's a more compact version

schro := {diff(psi(x), x, x)-(alpha*x^4+x^2-energy)*psi(x) = 0}

ic := {psi(3) = 0, (D(psi))(3) = 1}

schro1 := subs(energy = 3.30687, alpha = .1, schro)

soln1 := dsolve(schro1 union ic, {psi(x)}, type = numeric)

with(plots)

odeplot(soln1, [x, psi(x)], -3 .. 3)