Airy equation with ode solver

I am trying to solve Airy equation with matlab ode solvers (ode45, ode23, etc..)

Here is my function

function dPsidy = gpe(y,Psi)
    dPsidy = zeros(2,1);
    dPsidy(1) = Psi(2);
    dPsidy(2) = y.*Psi(1);
    

and my initial conditions

    yspan = linspace(5,-10,1000);
    Psi0 = [1e-4 1e-4];

After applying ode45 or other solvers I have solution that differs from built-in matlab airy function.

    [y,Psi] = ode45(@(y,Psi) gpe(y,Psi), yspan, Psi0);
    plot(yspan, Psi(:,1),yspan,airy(yspan))

How I can obtain the proper solution?

Answer 1

[Bjorn Gustavsson]

That is just a scaling-factor off. (in addition to some numerical issues around the zero-crossings.) This comes about because of your initial value. Try with the given values for airy(0) and its derivative there, i.e. 1/(3^(2/3)*gamma(2/3)) and -1/(3^(1/3)*gamma(1/3)). That you can integrate with ode45 in positive and negative direction - which will give you the proper comparison.

HTH

[anatolyb]

Well for the airy(0) it works, however I would like to extend this problem by adding nonlinearity term like so it is difficult to preddict what is happening at zero. I only know that for negative y the solution is close to but it didn't work.

Is there any general approach?

[Bjorn Gustavsson]

You misunderstand what I tried to explain. Your example comparison is wrong. That is due to your initial condition at x=5 is not consistent with airy(x), if you were to compare the ode45 solution with the proper initial condition you'll see that the ode45-solution is numerically fine:

x = [5-1/256,5,5+1/256];
fX = airy(x);
dfdX = gradient(fX,X);
psi0 = [fX(2),dfdX(2)] % compare with what you used
[y,Psi] = ode45(@(y,Psi) gpe(y,Psi), yspan, psi0);
plot(yspan, Psi(:,1),yspan,airy(yspan))

HTH