偏微分方程组有限差分法求解

问题遇到的现象和发生背景

我想问一下将一个偏微分方程组用有限差分法展开以后，如何将他转化为矩阵形式，最后代入matlab中求解？
目前偏微分方程我依据有限差分法已经展开了，但是由于他是一个四元二阶的方程组，我查阅了一些资料，均没有相关方面的解答，我看您发过类似方法的文章，还请您能不能提供一下大致的思路，亦或者是相关的文献参考，非常感谢！

图片中密度和速度是未知量，速度三个方向分别为uvw，这四个变量都是关于xyzt的函数

问题相关代码，请勿粘贴截图

运行结果及报错内容

我的解答思路和尝试过的方法

我想要达到的结果