def sturm_liouville_function(x, y, p, p_x, q, f, alpha=0, nonlinear_exp=2):
    """Second order Sturm-Liouville Function defining y'' for Lu=f.

    This form is used because it is expected for Scipy's solve_ivp method.

    Keyword arguments:
    x -- independent variable
    y -- dependent variable
    p -- p(x) parameter
    p_x -- derivative of p_x wrt x
    q -- q(x) parameter
    f -- forcing function f(x)
    alpha -- nonlinear parameter
    nonlinear_exp -- exponent of nonlinear term
    """
    y_x = y[1]
    y_xx = -1*(p_x/p)*y[1] + (q/p)*y[0] + (q/p)*alpha*y[0]**nonlinear_exp - f/p
    return [y_x, y_xx]