Numerical Integration after Numerical Solution of ODEs (MATLAB)

Question

My problem is similar to this one Numerical Integration, and I have already calculated the numerical solution of my ODEs using bvp4c, with the boundary conditions imposed at rmin=1e-5 (near r=0) and rmax=50. I do not have to keep the infinite interval, since even for $\frac{\lambda}{e^2}<<1$ the solution reaches the asymptotic behavior really fast.

So, I calculated my integral using trapz, but I would like to know if Matlab has a preciser way of doing this. I searched the methods user @drjrm3 mentioned in the question above, but I didn't understand what method I can implement when the integrand involves a combination of components of a vector which keeps the solution.

What I have so far is something like this:

f=trapz(xint,Sxint(3,:).^2. + 0.5*(1-Sxint(1,:).^2.).^2./xint.^2. + 0.5*xint.^2.*Sxint(4,:).^2. + ...
    Sxint(1,:).^2.*Sxint(2,:).^2. + 0.1*0.25*xint.^2.*(Sxint(2,:).^2. - 1).^2.)

Thanks in advance for any hint!

Answer 1

integral() and quadqk().

The exact method of integral() in Matlab changes over time. I haven't followed it for a while. A few years ago, I remember reading how it uses adaptive quadrature and were a slightly more advanced version of Guass-Kronrod. Now it seems Matlab doesn't talk about its integration method in the official documentation. Maybe they developed something good and proprietary. You can read the paper they linked the the documentation page and see if what they do is to your liking.