2 point BVP solver: how to compute errors

The way I computed the solution in the linked answer is the classical one: I just took a reference solution (assuming the code was correct, i.e. the numerical solution I found is the right one) with a small enough step size $h$, say $h=10^{-9}$. Then I computed the solution with smaller $h$s and for each one of those $h$ I computed the $|| e_h ||$ in a double logarithmic plot. I think in that case I used the supremum norm, because the choice of the parameters didn't produce boundary layers. As pointed out by Lutz in comments, this norm might not be the best one for that situations.

Be careful, this is telling you that the method is implemented correctly, not that your solution is the right one. If there are typos in your setup, you might be converging to another solution, not to the desired one.

I think you can use the Octave/MatLab code I provided in that answer to test the order of convergence. If you want to stick with Python, just define something like

m_range = [10**(i) for i in range(0,4)] #number of points in each grid

and loop over this. In the body, just compute the solution and check the error norm, appending this to a list, say error_list.

In the end, just do plt.loglog(m_range,error_list) or stuff like that. It's not important that your line fits the dots of the error, as in my answer (there I used a trick to shift it towards the error dots): the important point is that the error dots are parallel to the line with slope $-2$.