@sand15 Nice analysis. The wildly different orders of magnitudes of the various quantities makes me nervous, aside from the likely poorly posed problem.

@Alfred_F @janhardo Maple's isolve finds the general formula for the Pell's equation, so no extra procedure is required. (It expresses it in terms of an arbitrary integer parameter _Z1 rather than the equivalent recursion.) My point to @janhardo was that his procedure testing if integers on a grid satisfied this equation is not a realistic way since very large numbers can be needed even for equations with relatively small coefficients.

@janhardo The test-on-a-grid method fails fairly quickly - try to find a solution for x^2 - 151*y^2 = 1 aside from (1,0) and (-1,0).

It's not clear what you want to do "on line". If you just want help with your integration, please upload a worksheet showing the example and what you mean by unstable. Use the green up-arrow, choose your file, click upload and then one of insert link or insert contents.

@GunnerMunk If you put simplify in the definition of the velocity function, then it works.

This is a bug and I have submitted an SCR ("bug report") to Maple.

Download Error_2.mw

@MichaelVio Here is the other approach. I don't follow what you are trying to do, but the values you are using don't seem reasonable. Note that the approximations made in deriving the Planck law assume the volume is much much larger than the photon wavelengths and there are many many photon energy states in the frequency range; not sure they apply here.

phot3.mw

@GunnerMunk Please upload.

@C_R I think the new rule based system is likely an extra layer on top of the existing system to give more flexibility, but didn't change the existing lower-level system.

@salim-barzani It works after you remove "=0".

sys-Dr-D.mw

@salim-barzani As usual I do not set the ODE =0.
If you put the coeffs in a list rather than a set, they are in 1:1 correspondence with the monomials.

system1.mw

@MichaelVio The worksheet uses dsolve to solve your equation, but it has a singularity and won't work for the values you use.

phot2.mw

It explains how to find k*T, but doesn't work here because of the strange function.

In terms of the photon gas, it seems what you are doing is backwards. In the photon gas we have a formula for diff(E(nu),nu) which is a function g(nu) of nu (or omega) Eq 63.5, which is a peaked function (Fig. 7). So to find E(nu) you would integrate the right-hand side. If you want h*nu on the right hand side to be a function h(nu)*nu or E(nu) because you are doing some different physics, then the differential equation would just be diff(E(nu),nu) = g(nu,E(nu)).

If one integrates a peaked function, then the energy will be an increasing function of nu, and you will not find two values of nu with the same energy.

Instead you have written E(nu) = diff(g(nu,E(nu)), nu) for your similar function g, which seems backwards to me.

Edit: I am thinking of E(nu) as the energy of the photons in the gas with frequencies up to nu, but maybe you want some other interpretation. I am not following what you are trying to do.

@MichaelVio I cannot see the Landau&Lifsch.pdf, probably because the system doesn't like the & in the filename. Please upload again.

For 1, you have both he(nu) and h(nu), but I assume you mean to solve the de that led to it, assuming he(nu)=h(nu).

For 2, I am assuming that when you say "fit" you mean find the single value of k*T that makes those two conditions true.

@MichaelVio You said:

"Please advise on solving the Planck.mw and plot the solution for h(nu) with E=ν⸱h(nu). Thus, for nu2 =9.733521364*10^16 =>h(nu2)= 3.348222989*10^(-17) eV*s and the same energy of 3.259eV for nu1 = 7.889275211*10^14 =>h(nu2)=4.135667697*10^(-15) eV*s  With the partular case the plot 

with the supposition that k = 1/2 eV/K, thus kT=>1.380649*10^(-23)*297/(1.60217*10^(-19));
=0.02559, so kT should be the energy close to 0.5 and could be assimilated to an energy E0"

This is my understanding of this: You have a differential equation in h(nu) for which you want the energies at both nu1 = 7.889275211*10^14 s^(-1) and nu2 =9.733521364*10^16 s^(-1) to be the same energy E1 = 3.259eV. You want to find the value of k*T that makes this true. This can be done, though there might be numerical issues.

Since k is just Boltzmann's constant, it seems you want to find the temperature for which this is true? But then I don't know why you specify T = 297 K.

If I open to my default blank worksheet and then close the tab so that there is no open worksheet, then I can see the recent documents list but trying to open the most recent document fails (nothing happens), but I can open other recent documents further down the list.

I had noticed this before in the context of working on the most recent document, then closing its tab to leave no open worksheets, then I cannot reopen it from the recent documents list.

@MichaelVio So I agree with your dimensional analysis. At some point you set k:=1/2. But k is Boltzmann's constant. I think you want to change Boltzmann's constant to be k__B, and then k is something else. But I don't understand what you want to do.