用python数值求解一维偏微分方程组（pypde lib或任何其他）

Question

我有以下偏微分方程组：

在哪里 c0-常数， r-独立空间变量， t-时间变量， f（r，t）-第一个未知函数， f_t（r，t）-第二个未知函数（实际上它只是表示f对t的一阶导数：f_t（r，t））， f_r（r，t）-f对r的一阶导数， f_rr（r，t）-f对r的二阶导数

必须在变量范围r0..r1，t0..t1（r0=5*10^-4，r1=1，t0=0，t1=5*10^-5）内找到未知的f（r，t）和f\u t（r，t）

还有以下初始条件：

f(r, t0) = 0
f_t(r, t0) = 0

以及边界条件：

f_t(r0, t) = ...<a function which is defined by set of points>...
# e.g.: scipy.interpolate.interp1d(array_of_t_values, array_of_F_values, kind='linear')

如何定义所有这些，并用python中的pypde（或任何其他合适的库）进行数值求解？我还需要画出f（r，t）和f_t（r，t）的结果曲面

下面我试图修改现有示例（py pde），但我不知道如何定义导数：

from pde import FieldCollection, PDEBase, UnitGrid

class MyPDE(PDEBase):
    def __init__(self, c0=1500):
        self.c0 = c0

    def evolution_rate(self, state, t=0):
        f, f_t = state

        f = f_t
        f_t = (f.diff(r, 2) + 1/r * f.diff(r)) / self.c0

        return FieldCollection([f, f_t])


grid = UnitGrid([32, 32])
state = FieldCollection.scalar_random_uniform(2, grid)

eq = MyPDE()
result = eq.solve(state, t_range=100, dt=0.01)
result.plot()