Numerically solving a partial differential equation in python with Runge Kutta 4

Question

I'm supposed to solve the following partial differential equation in python using Runge-Kutta 4 method in time.

$$ \frac{\partial}{\partial t}v(y,t)=Lv(t,y) $$

where $L$ is the following linear operator

$$ L=\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)^{-1}\left\{-i\alpha(1-y^2)\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)-2i\alpha+\frac{1}{R}\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)^2\right\} $$

The spatial dimension is discretized along the $N = 101$ grid points from $y = -1$ to $y = 1$.

The boundary conditions are v[0] = 0 and v[N-1] = 0 and

$$ \left.\frac{\partial v}{\partial y}\right|_{y=-1}=\left.\frac{\partial v}{\partial y}\right|_{y=1}=0 $$

I have to discretize the second- and fourth-derivative operators using central finite difference and create a function for them. Using those functions, I can then create the $L$ operator.

Finally, using the initial condition $v(y) = 0.02(1+\cos(\pi y))$ time step $dt= 0.01$ and parameter values $R=500$ and $\alpha=0.3$, I have to find $v$ after 10s.

I have written the following code where the $L$ operator is put in a matrix form. The first few iterations seem to work giving somewhat reasonable results but then the values explode and I get an error. I'm not sure where the problem is, I might be misunderstanding how to implement the problem.

Another point I'm unsure about is whether the inverse of the $\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)$ operator can be understood as the inverse of the matrix representation of that operator or if that's not the case.

Can anybody give me some pointers as to how to solve this?

import numpy 
import scipy
from scipy import linalg
from matplotlib import pyplot
%matplotlib inline

N=101
L=1.0

#creating the grid points
y=numpy.linspace(-L+0.04,L-0.04,N-4)
dy=(2*L)/(N-1)

v=numpy.empty(N-4)

# Operator of the second derivative acting on v with respect to y 
def D2_v(N,dy):

    # Setup the diagonal of the operator.
    D = numpy.diag((-2.0/ (dy**2)) * numpy.ones(N-4))

    # Setup the upper diagonal of the operator.
    U = numpy.diag((1.0/ (dy**2)) * numpy.ones(N - 5), k=1)
    # Setup the lower diagonal of the operator.
    L = numpy.diag((1.0/ (dy**2)) * numpy.ones(N - 5), k=-1)
    # Assemble the operator.
    D2_v = D + U + L
    # Setup the Neumann condition 
    D[-1, -1] = (-7.0/4.0) / (dy**2)
    D[0, 0] = (-7.0/4.0) / (dy**2)

    return D2_v

# Operator of the fourth derivative acting on v with respect to y 
def D4_v(N,dy):

    # Setup the diagonal of the operator.
    D = numpy.diag((6.0/ (dy**4)) * numpy.ones(N-4))

    # Setup the upper diagonal of the operator.
    U = numpy.diag((-4.0/ (dy**4)) * numpy.ones(N - 5), k=1)
    U2 = numpy.diag((1.0/ (dy**4)) * numpy.ones(N - 6), k=2)

    # Setup the lower diagonal of the operator.
    L = numpy.diag((-4.0/ (dy**4)) * numpy.ones(N - 5), k=-1)
    L2 = numpy.diag((1.0/ (dy**4)) * numpy.ones(N - 6), k=-2)

    # Assemble the operator.
    D4 = D + U + L + U2 + L2

    # Setup the Neumann condition 
    D4[-1, -1] = (5.0) / (dy**4)
    D4[0, 0] = (5.0) / (dy**4)
    D4[1, 0] = (-7.0/4.0) / (dy**4)
    D4[-2, -1] = (-7.0/4.0) / (dy**4)

    return D4

D2=D2_v(N,dy)
D4=D4_v(N,dy)   

# Operator L
def L_v(N,y,dy,R,alpha):

    Ma = alpha*numpy.identity(N-4)
    My = numpy.diag(y)
    I = numpy.identity(N-4)
    Lv = numpy.linalg.inv(D2-Ma**2)*(-1j*Ma*(I-My**2)*(D2-Ma**2)+1j* 
     (-2.0)*Ma+(1.0/R)*(D4-D2*(Ma**2)-(Ma**2)*D2+(Ma**4)))

    return Lv    

#parameters
alpha=0.3
R=500.0
dt=0.01

L1=L_v(N,y,dy,R,alpha)

#initial value for v
for i in range(N-4):    
    v[i]=0.02*(1.0+numpy.cos(numpy.pi*y[i]))

#Runge Kutta implementation for 1000 time steps
for m in range(1000):

    k1 = dt*numpy.matmul(L1,v)
    k2 = dt*numpy.matmul(L1,v + 0.5*k1)
    k3 = dt*numpy.matmul(L1,v + 0.5*k2)
    k4 = dt*numpy.matmul(L1,v + k3)

    v = v + (1.0/6.0)*(k1+2*k2+2*k3+k4)