Converting symbolic ODE solution to function ... Warnings

I am trying to find symbolic solution of specific non-linear ODE and then convert it to the matlab function via flowing code:

syms y(x) f(x) 
syms x x0 y0 a b c d real
eqn = diff(y,x) == f/(a+b*y+c*y^2+d*y^3);
cond = y(x0)==y0;
S(x) = dsolve(eqn,cond);
sS(x) =simplify(S)
matlabFunction(sS,"File","myfunc.m");

where

a,b,c,d are scalar real constants

f(x) is known external real function

x0,y0 define initial condition of ODE problem

During matlab code conversion I always get a several Warnings:

Warning: symbolic:generate:FiniteSetsRootNotSupported∦'root' without index not supported. Using index '1' instead.∦
Warning: symbolic:generate:FunctionNotVerifiedToBeValid∦Function 'f' not verified to be a valid MATLAB function.∦
Warning: symbolic:generate:FiniteSetsRootNotSupported∦'root' without index not supported. Using index '1' instead.∦
Warning: symbolic:generate:FunctionNotVerifiedToBeValid∦Function 'f' not verified to be a valid MATLAB function.∦
Warning: symbolic:generate:FunctionNotVerifiedToBeValid∦Function 'f' not verified to be a valid MATLAB function.∦
Warning: symbolic:generate:FiniteSetsNotSupported∦Finite sets ('DOM_SET') not supported. Using element '(t4 + a*y0)/a' instead.∦
Warning: symbolic:generate:FiniteSetsNotSupported∦Finite sets ('DOM_SET') not supported. Using element 't14*(t6 + t20)' instead.∦
Warning: symbolic:generate:FiniteSetsNotSupported∦Finite sets ('DOM_SET') not supported. Using element '-t14*(a - (t7 + t13 + t16 + t19)^(1/2))' instead.∦
Warning: symbolic:generate:FiniteSetsNotSupported∦Finite sets ('DOM_SET') not supported. Using element '-t14*(a + (t7 + t13 + t16 + t19)^(1/2))' instead.∦

and the following matlab code:

function sS = BURh103(x,a,b,c,d,x0,y0)
%BURh103
%    sS = BURh103(X,A,B,C,D,X0,Y0)
%    This function was generated by the Symbolic Math Toolbox version 24.2.
%    03-Oct-2024 15:16:31
t0 = roots([d.*-3.0,c.*-4.0,b.*-6.0,a.*-1.2e+1,a.*y0.*1.2e+1+b.*y0.^2.*6.0+c.*y0.^3.*4.0+d.*y0.^4.*3.0+integral(@(u)f(u),x0,x).*1.2e+1]);
t2 = t0(1);
t0 = roots([c.*-4.0,b.*-6.0,a.*-1.2e+1,a.*y0.*1.2e+1+b.*y0.^2.*6.0+c.*y0.^3.*4.0+integral(@(u)f(u),x0,x).*1.2e+1]);
t3 = t0(1);
t4 = integral(@(u)f(u),x0,x);
t5 = b.*t4;
t6 = b.*y0;
t7 = a.^2;
t9 = (b ~= 0.0);
t10 = (c == 0.0);
t11 = (d == 0.0);
t14 = 1.0./b;
t15 = sqrt(2.0);
t13 = t5.*2.0;
t16 = a.*t6.*2.0;
t17 = a+t6;
t18 = sqrt(t5);
t19 = t6.^2;
t20 = t15.*t18;
if ~all(cellfun(@isscalar,{b,c,d,t10,t11,t17,t9}))
    error(message('symbolic:sym:matlabFunction:ConditionsMustBeScalar'));
end
if (d ~= 0.0)
    sS = t2;
elseif (t11 & (c ~= 0.0))
    sS = t3;
elseif ((t10 & t11) & (b == 0.0))
    sS = (t4+a.*y0)./a;
elseif (((t9 & t10) & t11) & (t17 == 0.0))
    sS = t14.*(t6+t20);
elseif (((t9 & t10) & t11) & (0.0 < t17))
    sS = -t14.*(a-sqrt(t7+t13+t16+t19));
elseif (((t9 & t10) & t11) & (t17 < 0.0))
    sS = -t14.*(a+sqrt(t7+t13+t16+t19));
else
    sS = NaN;
end
end

I did not understand the following warning messages:

Warning: symbolic:generate:FiniteSetsRootNotSupported∦'root' without index not supported. Using index '1' instead.∦
Warning: symbolic:generate:FiniteSetsRootNotSupported∦'root' without index not supported. Using index '1' instead.∦
Warning: symbolic:generate:FiniteSetsNotSupported∦Finite sets ('DOM_SET') not supported. Using element '(t4 + a*y0)/a' instead.∦
Warning: symbolic:generate:FiniteSetsNotSupported∦Finite sets ('DOM_SET') not supported. Using element 't14*(t6 + t20)' instead.∦
Warning: symbolic:generate:FiniteSetsNotSupported∦Finite sets ('DOM_SET') not supported. Using element '-t14*(a - (t7 + t13 + t16 + t19)^(1/2))' instead.∦
Warning: symbolic:generate:FiniteSetsNotSupported∦Finite sets ('DOM_SET') not supported. Using element '-t14*(a + (t7 + t13 + t16 + t19)^(1/2))' instead.∦

What exactly these messages are trying to say?

Moreover, the matlab code is far from optimal or effective code. The auxilliary variable txx are precomputed, but the same integral(@(u)f(u),x0,x) is computed three times?!

[Torsten]

Your differential equation simply has no explicit analytical solution that could be evaluated in a MATLAB function. The main problem is the general dependence of f on x. But in either case, the solution will be implicit rather than explicit in y, i.e. something like F(x) - G(y) = 0 where the inverse of G is difficult to determine (G is a polynomial of degree 4 in y).

[Michal]

Yes no explicit solution, but possible numerical solution by roots function for each separate x.

Answer 1

[Torsten]

The solution is implicitly given by the equation

(a*y+b*y^2/2+c*y^3/3+d*y^4/4)-(a*y0+b*y0^2/2+c*y0^3/3+d*y0^4/4) - integral_{x'=x0}^{x'=x} f(x') dx' = 0

So given a specific value for x, you could in principal evaluate the integral over f with limits x0 and x and solve the resulting polynomial in y to get four different solutions (from which you had to select the correct one).

But my advice is to solve the equation numerically using ode45, e.g.