# Help with code for FDM Couette flow of two fluids on python

I'm trying to solve numerically with a python code a very simple BVP that represent the velocity field of two fluids between two plates where the driving force of the momentum transport is only for the movement either upper plate or bottom plate. $$V$$ is the velocity (not in module) of bottom and $$U$$ is the velocity of upper plate. $$m1$$ and $$m2$$ are the viscosity of fluids. Here the equations:

$$\frac{d^2 v(y)} {dy^2} = 0$$

with:

$$v(0)=V$$

$$v(\varepsilon H)=u(\varepsilon H)$$

And:

$$\frac{d^2 u(y)} {dy^2} = 0$$

with:

$$v'(\varepsilon H) m1=m2 u'(\varepsilon H)$$

$$u(H)=U$$

To solve the problem I tried to use the centered scheme of FDM at #2 order of accuracy (there isn't for first order acc.), and for the Neumann condition I tried to use for $$dv/dy$$ the backward scheme and for %du/dy\$ the forward scheme. The problem is the interface, I cant' find a numerical solution because the problem doesn't converge.

This is my python code, i hope that there are an error:

import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

# Parameters
e1 = 1 #viscosities
e2 = 5
H = 2
h = 0.001
U = 14
V = 0
ep = 0.3 # ratio between thickness of fluid 1 and total thickness

#nodes
N = int(H/h+1)
a = e1/e2
M = int(ep*H/h+1) #interface node

def eq_field(z):

v = z[0:M+1] #da 0 a M
u = z[M+1:] #da M+1 fino alla fine
f = np.zeros(N+1)

f[0]= -V + v[0]

for i in range(1,M):
f[i]=v[i-1]-2*v[i]+v[i+1]

f[M] = v[-1]-u[0]
f[M+1] = v[-1]-e2/e1*(u[1]-u[0])-v[-2]
f[M+2] = e2/e1*(u[1]-u[0])+v[-2]-2*u[1]+u[2]

for i in range(1,N-M-3):
f[i+M]=u[i-1]-2*u[i]+u[i+1]

f[-1]=u[-1]-U

return f

p1 = np.linspace(0,U/2,M+1)
p2 = np.linspace(U/2,U+1,N-M)
p = np.concatenate((p1,p2))

field = fsolve(eq_field,p)
print(field[int((ep*H)/h+1)])
print(field)
x = ep*H*np.ones(len(field))
y = np.linspace(0,H,len(field))
plt.plot(field,x)
plt.plot(field,y)
plt.show()'


I appreciate every tip.