def eulers_richardson_method(f, dx, y, yp, range, return_yp = False):
    """ The Eulers Richardson method for solving a differential equation of second
    order. Works by taking the Euler method, but uses the average values instead.
    This produces a better approximation to the solution faster (in theory) than
    just the plain Eulers method.
    
    :param f: The input math function derivative whom to approximate its solution
    :type f: lambda
    :param dx: The step size to use. Smaller is better.
    :type dx: float
    :param y: The initial value of y given by initial condition.
    :type y: float, int
    :param yp: The initial value of y' given by initial condition.
    :type yp: float, int
    :param range: A list which specifies the beginning and the ending of our domain.
    :type range: list
    :return: Returns a tuple for the x coordinates corresponding to a set of y coordinates,
    which approximate the solution to f.
    :rtype: list
    """
    x       = min(range)
    y_space = [y]
    yp_space = [yp]
    x_space = [x]
    
    while x<=max(range):
        yp_mid  = yp + 1/2*f(x,y,yp)*dx
        y_mid   = y  + 1/2*yp*dx
        ypp_mid = f(1/2*x*dx, y_mid, yp_mid)
        yp      += ypp_mid*dx
        y       += yp_mid*dx
        
        x       += dx
        x_space.append(x)
        y_space.append(y)
        yp_space.append(yp)
    if (return_yp):
        return (x_space, y_space, yp_space)
    return (x_space, y_space)