Python solve_bvp fourth order differential equation

Question

I'm trying to apply scipy's solve_bvp to the following problem

T''''(z) = -k^4 * T(z)

With boundary conditions on a domain of size l and some constant A:

T(0) = T''(0) = T'''(l) = 0

T'(l) = A

So far, I've reduced the fourth order equation to a 1st order system and written the following function:

def fun1(t, y):
    y0 = y[1]
    y1 = y[2]
    y2 = y[3]
    y3 = -k**4 * y0
    ret = np.vstack((y0, y1, y2, y3))
    return ret

Then, I've established my boundary conditions, trying to follow the documentation (which I don't really understand...)

def bc(ua, ub):
    # 0th, 1st, 2nd and 3rd derivative BCs
    return [ua[0], ub[1]-A, ua[2], ub[3]]

Then I set up my initial guesses

A, l = 10, 3
x_init = [0, l]
y_init = [[0, 0], [0, A], [0, 0], [0, 0]]

when I run solve_bvp(fun, bc, x, y), however, I get the wrong solution. I don't know quite why. The solver converges, but it doesn't look like how I expect it to.

Can someone please explain what the bc function is supposed to return for Von Neumann boundary conditions? I'm really struggling to wrap my head around the documentation...

Answer 1

In computing y3, you need to actually use y[0], not y0=y[1].

To avoid such misunderstanding I would write

def fun1(t, y):
    dy0 = y[1]
    dy1 = y[2]
    dy2 = y[3]
    dy3 = -k**4 * y[0]
    return np.vstack((dy0, dy1, dy2, dy3))