def shooting_dirichlet(f, interval, bcs, N):

    """
    Solve a BVP using shooting

    Parameters
    ----------

    f: function
        RHS function defining BVP: y'' = f(x, y, y')
    interval: list of two floats
        [a, b] where BVP defined on x \in [a, b]
    bcs: list of two floats
        Defines Dirichlet boundary values [y(a), y(b)]
    N: single integer
        Number of interior points in the grid, used for output only

    Returns
    -------

    x: array of floats
        Location of grid points
    y: array of floats
        Solution at grid points
    """

    s = np.linspace(interval[0], interval[1], N+2)
    y = np.zeros_like(s)
    y[0] = bcs[0]

    def rhs_shooting(s, z):
        """
        RHS function for shooting (z = [y, y']; solve z' = [z[1], f(...)])

        Parameters
        ----------

        x: float
        Location
        z: list of floats
        [y, y']

        Returns
        -------

        dzdx: list of floats
        [y', y''] = [z[1], f(x, z[0], z[1])]
        """

        v = z[1]
        theta = z[0]
        dtheta = v
        dvbyds = s * fg * np.cos(theta) + s * fx * np.sin(theta)
        #y, dy = z
        return np.array([dtheta, dvbyds])

    def residuals_shooting(guess):
        """
        Solve IVP given y'(a): compute r = y(b) - bcs[1]

        Parameters
        ----------

        guess: float
        Guess for derivative of solution at x=a

        Returns
        -------

        residual: float
        Values of residual r at x=b
        """
        sol_shoot = solve_ivp(rhs_shooting, interval, [bcs[0], guess])
        return sol_shoot.y[1, -1] - bcs[1]


    sol = scipy.optimize.root_scalar(residuals_shooting, method='brentq', bracket=[-10, 10])
    guess = sol.root
    sol_shoot = solve_ivp(rhs_shooting, interval, [bcs[0], guess], t_eval=s)
    y_all = sol_shoot.y
    return s, y_all[0, :]

        #Solve IVP with theta
dxds = np.cos(y_all)
dzds = np.sin(y_all)
x0 = R*np.cos(theta0)
z0 = R*np.sin(theta0)

x = solve_ivp(dxds, interval, x0)
z = solve_ivp(dzds, interval, z0)