# Solve a system of linked differential equations with solve_bvp

I have to solve the following problem :

I wrote the second hand :

def dN_dthyb(t,N,param):

alpha=param[:-3]
r0=param[-1]
m=len(alpha)+1
dN = np.zeros(m)
dN[0] = r0-alpha[0]*N[0]
for i in range(1, m-1):
dN[i] = alpha[i-1]*N[i-1]-alpha[i]*N[i]
dN[m-1] = alpha[m-2]*N[m-2]
return dN


Where param is an array with the values of alpha and the value of r(t) (here r is constant) I already used this function in order to solve the problem as a Cauchy problem with ode and now I want to have the boundary coundition such as N(0) = N(T), and r is a periodic function with amplitude equal to r0.

First, in order to use solve_bvp, I rewrore the function :

def dN_dthyb(t,N,param):
alpha=param[:-3]
beta=param[-(2+nbr_maturation):-3]
r0=25
h1=param[-2]
h2=param[-1]
m=len(alpha)+1
n=len(beta)+1
dN = np.zeros((m+n,len(t)))
dN[0,:]=r(t,r0)-alpha[0]*N[0,:]
for i in range(1, m-1):
dN[i,:] = alpha[i-1]*N[i-1,:]-alpha[i]*N[i,:]
dN[m-1,:] = alpha[m-2]*N[m-2,:]
dN[m,:] = -beta[0]*N[m,:]
for i in range(1, n-1):
dN[m+i,:]=beta[i-1]*N[m+i-1,:]-beta[i]*N[m+i,:]
dN[-1]=beta[-1]*N[-2,:]
return dN


Is it correct to see N as a matrix where each row, for example row i represents N_i(t) ?

Secondly, I don't understand how to set the boundary condition. I don't understand the nature of the return value. What does this value correspond to?

• from scipy.integrate import solve_bvp, note the "from". Is there a reason to suspect that periodic solutions exist? You would need strong forcing terms $p(t)$, $r(t)$ to overcome the mostly exponential behavior from the remaining terms Aug 4 at 14:47