odepkg would be the logical candidate, but AFAIK it hasn't been maintained for several years

bvp4c (along with bvpinit) is a core Matlab function currently not in Octave. So it belongs in core Octave.

Fixed a typo in the finite element method which would cause a crash if the number of equations does not equal 2.

(file #50905)

I have implemented bvp4c using a finite difference method which is different from the Lobatto IIIA methods used in Matlab.

The Theoretical Basis of the Method

In my method to obtain an estimate of the derivative, I use the nearest five points to calculate the coefficients of a fourth-order polynomial passing through the potentially unevenly spaced points. This part of the method in the subfunction finitedifferencecoef1d follows the normal way of deriving finite-difference rules with more than two points. It matches coefficients of the fit with the Taylor series expansion about the point we want to take the derivative.

I use the inbuilt solver fsolve to obtain the values of the function with a residue calculated from my finite difference derivative minus the user-supplied function. Boundary conditions are appended.

I estimate the error by comparing my functions derivative with the user-supplied derivative at the midpoint between two adjacent points. If the difference is greater than the tolerance I add that point and resolve until all midpoints are within the required tolerance.

Why Not to use Matlab's Method

The Lobatto IIIA method used by Matlab is an implicit Runge–Kutta method. Runge–Kutta methods estimate the value of the function using recursive substitution of weighted derivatives (different locations) into the derivative function. This recursive substitution requires loops within the residue function which may be called many times during the solution process. This inability to vectorize the residue will make any pure Octave implementation very slow. I found an implementation that I think was once part of Octave's ode package that uses the Lobatto IIIA method. This implementation is unbearably slow.

Comparison of the Methods

The main advantage of my method is the implementation is almost completely vectorized. The residual function is completely vectorized. Profiling the runtime shows almost all of the time is spent in the solver with \ the most time-consuming operation.

My method has no subpoints between the points returned by the method. This results in a larger system of equations needing to be solved.

As the method uses the nearest five points to calculate the function and derivative, the function and the derivative are not continuous.

Current State of the Code

The code can solve odes efficiently. Two demos included.

It is missing the following:

• Ability to apply constraints in the middle of the domain
• Solve for unknown parameters
• Use analytical derivatives

My method is bvp4cfd.m it uses bvpinit.m to prepare the initial guess. The method I found which I think was part of the ode package is bvp4cold.m. (bvp4cold.m is not my code and I am not suggesting to include it. It highlights the difficulties of using the Runge–Kutta method.)

I have not prepared the documentation or followed the style guide I will fix these issues after I finalize the code.

Information I Would Like Before Finishing the Code

Is it acceptable to use a different method to what Matlab uses to solve the same problem?

I'm planning to write a finite element method for bvp4c which will allow a continuous function and first derivative. It contains subpoints so the system of equations the solver needs to solve should be smaller. I think vectorized code should be possible. If I can get the finite element method to a comparable speed to this implementation would there be a preference for finite difference or finite element method?

Let me know if any of my input or output is different from Matlab.