kill(all);
load('contrib_ode);

batch("transform.max");

batch("determine_coeff.max");

batch("ode2sing_2020_Jan_23.max");

## Parabolic cylinder functions

ode2sing can compute solutions in terms of parabolic cylinder functions. The simplification of the result should be improved, e.g. constant factors should be eliminated.

p1: 0;
p0: -x^2+x+1;
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

The Hermite differential equation:

kill(lambda);
lambda: 7;
p1: -2*x;
p0: 2*lambda;
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

Solution by a composition of Kummer function with a power function:

kill(a);
p1: 0;
p0: -(x-8)^(2*a-2)+(x-8)^(a-2);
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

## Mathieu equation and spheroidal wave equation

ode2sing can find solutions which are spheroidal wave functions and also checks whether they are Mathieu functions.

a: 3;
q: 5;
p1: -x/(1-x^2);
p0: (a+q*(2-4*x^2))/(1-x^2);
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

b: 2;
c: 6;
q: 8;
p1: -2*(b+1)*x/(1-x^2);
p0: (c-4*q*x^2)/(1-x^2);
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

## Bessel equation and Kummer equation

ode2sing can find solutions which are the composition of Bessel or Kummer functions with a power of a Mobius transformation.

p1: 0;
p0: (4*(x-1)^4-15*(x-3)^4)/(4*(x-1)^2*(x-3)^6);
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

The Bessel equation:

a: 5;
p1: 1/x;
p0: 1-a^2/x^2;
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

The Kummer equation: ode2sing can find solutions which are Laguerre functions (a is a negative integer)

a: -5;
b: 1;
p1: b/(x-2)-1;
p0: -a/(x-2);
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

A Bessel function applied to a fractional power:

p1: 0;
p0: (36*x^(2/3)+5)/(36*x^2);
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

The following check does not return 0. I am not sure whether this indicates an error in my program.

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

An example in which ode2sing finds the solution with less computation time than contrib_ode

p1: 1/x;
p0: 16*x^6-16/x^2;
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

contrib_ode(eq,y,x);

## Hypergeometric equation

ode2sing can find solutions which are Legendre functions.

mu: 2;
nu: 3;
p1: 2*x/(x^2-1);
p0: (-nu*(nu+1)*x^2+nu*(nu+1)-mu^2)/(x^2-1)^2;
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

ode_check(eq,legend[1]);

contrib_ode(eq,y,x);

In the following example contrib_ode gives a more explicit solution than ode2sing since ode2sing does not use some special cases for hypergeometric functions.

a: 2;
b: 5;
c: 2;
p1: (c-(a+b+1)*x)/(x*(1-x));
p0: -a*b/(x*(1-x));
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

A hypergeometric equation with complex parameters:

p1: 1/(1+x^2);
p0: 1/(1+x^2);
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

An example of an equation with 4 regular singular points:

p1: 0;
p0: (x^2+1)/(4*(x+1)^2*x^2*(x-1)^2);
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p1,p0,y,x);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

## Transformations

ode2sing can find the solutions to the Mathieu equation if the user suggests a suitable coordinate transformation and its inverse:

p1: 0;
p0: a-2*q*cos(2*x);
eq: 'diff(y,x,2)+p1*'diff(y,x)+p0*y=0;
ode2sing(p0,p1,y,x);

phi: cos(x);
phi_inv: acos(x);
ode2sing_transform(p1,p0,y,x,phi,phi_inv);

ode_check(eq,%[1]);

contrib_ode(eq,y,x);

The following equation describes the theta-part of the spherical harmonic functions. ode2sing can find the solutions if the user suggests a suitable coordinate transformation and its inverse.

kill(ell);
kill(m);

p1: cos(theta)/sin(theta);
p0: ell*(ell+1)-m^2/sin(theta)^2;
eq: 'diff(y,theta,2)+p1*'diff(y,theta)+p0*y=0;
ode2sing(p1,p0,y,theta);

phi: cos(theta);
phi_inv: acos(theta);
ode2sing_transform(p1,p0,y,theta,phi,phi_inv);

The solution in terms of Legendre functions is not strictly correct. One has to replace the absolute value of m by m.

ode_check(eq,%[1]);

contrib_ode(eq,y,theta);