3.19 How to find the coefficients of any ode, linear or not and any order?

This function below finds the ode coefficients for any ode, linear or not. It takes as input the ode and returns back two items. The first is list of coefficients, and the second is the RHS (i.e. \(f(x)\) ).

So give an ode such as \( 3 y'' + 5 y' =\sin x\) it will return [0,5,3],sin(x) where the order of coefficients is \(y,y',y'',\dots \). If terms is missing, zero is placed there.

#version Dec 3, 2024, 3PM. by Nasser M. Abbasi 
get_ode_coefficients_any_order:=proc(ode_in::`=`,func::function(name),$)::list,anything;  #a,f 
 
#find coefficients. For example, for 5*y''''+y'''+y'+3*y=sin(x); it will return 
#   [3,1,0,1,5],sin(x) 
#Where the list is in order y,y',y'',y''',y'''' 
 
local y:=op(0,func); 
local x:=op(1,func); 
local ode:=ode_in,i,j,tmp; 
local a::list; 
local c:=0; 
local f:=0; 
local ode_order::integer; 
local found_diff::truefalse; 
 
    ode:= expand(lhs(ode))-expand(rhs(ode)); 
    ode_order:= PDEtools:-difforder(ode,x); 
    if ode_order<1 then 
        error "ODE expected in get_ode_coefficients_any_order, found ", ode; 
    fi; 
 
    a:=[0$ode_order+1]; #for [y,y',y'',...,y(n)], that is why added 1 
 
    if ode::`+` then 
        for i from 1 to nops(ode) do 
            tmp := op(i,ode); 
            found_diff:=false; 
            for j from ode_order to 1 by -1 do 
                if has(tmp, diff(y(x),x$j) ) then 
                    a[j+1] := a[j+1] + tmp/ diff(y(x),x$j); 
                    found_diff:=true; 
                    break; 
                fi; 
            od; 
            if not found_diff then 
               if has(tmp,y(x)) then 
                    a[1] := a[1] + tmp/y(x); 
                else 
                    f := f+tmp; 
                fi; 
            fi; 
        od; 
    else 
        found_diff:=false; 
        for j from ode_order to 1 by -1 do 
            if has(ode, diff(y(x),x$j) ) then 
               a[j+1] := ode/ diff(y(x),x$j); 
               found_diff:=true; 
               break; 
            fi; 
        od; 
        if not found_diff then 
            if has(ode,y(x)) then 
                a[1] := ode/y(x); 
            else 
                f :=ode; 
            fi; 
        fi; 
    fi; 
 
    return a,(-f); 
 
end proc:
 

Examples

ode:=sin(x)+y(x)*diff(y(x),x$5)+diff(y(x),x)-1/x+99*diff(y(x),x$3)=0; 
get_ode_coefficients_any_order(ode,y(x)) 
                                                1 
              [0, 1, 0, 99, 0, y(x)], -sin(x) + - 
                                                x 
 
ode:=diff(y(x),x$3)=0: 
get_ode_coefficients_any_order(ode,y(x)) 
                        [0, 0, 0, 1], 0 
 
ode:=diff(y(x),x$3)=sin(y(x)): 
get_ode_coefficients_any_order(ode,y(x)) 
 
                   [  sin(y(x))         ] 
                   [- ---------, 0, 0, 1], 0 
                   [    y(x)            ] 
 
 
ode:=sin(x)+1/x*diff(y(x),x$3)+y(x)=exp(x)+9; 
get_ode_coefficients_any_order(ode,y(x)) 
 
               [         1] 
               [1, 0, 0, -], -sin(x) + exp(x) + 9 
               [         x]