# Solving for Coefficients in a System of Linear Equations Using Mathematica

I'm currently using Matlab to solve a two-dimensional PDE in a rectangular domain using finite differences. Without going into too much detail (although I'd be happy to if necessary), I have the following discretized equation: $$u_{j,k}^{n+1}-\Delta t\ L_{j,k}^{n+1}=u_{j,k}^n \qquad\text{for } j=1,...,N-1\ \text{ and } k=0,...,M$$
where $$L_{j,k}^{n+1}$$ represents the discretized scheme for the PDE and involves a seven-point stencil (three neighboring points in $$x$$ and five neighboring points in $$z$$, both including the centered point). I used the relevant boundary conditions in $$z$$ to determine equations for the points just outside my domain (i.e. at $$u_{j,-2},u_{j,-1},u_{j,M+1},u_{j,M+2}$$) in terms of points inside my domain. The boundary conditions in $$x$$ are Dirichlet conditions which is why in the above equation $$j$$ only ranges from $$1$$ to $$N-1$$.

I have been, by hand, attempting to arrange these equations into a system of linear equations so Matlab can solve for all variables at once by mapping the points in 2D space to a 1D vector by: $$l =(j-1)(M+1)+k+1$$ I was then hoping to arrange the equations into the form: $$Q_1\vec{u}^{n+1}=Q_0\vec{u}^n+\vec{b}$$ Since I'm using a forward Euler scheme, $$Q_0=I$$ and $$\vec{b}$$ is only non-zero near the boundaries in $$z$$ (due to the specific BCs applied).
However, my code has some bugs and I'm assuming it's from simple mistakes made during this process. Is it possible to have Mathematica do this for me? I know this is a quite involved question so any example or other resource you could link would be greatly appreciated. If it helps, I can include the equations for the points outside the domain. I am interested in knowing how to do something like this on a computer to limit mistakes in all my studies going forward though, so a general example would also be just as good.

EDIT: The PDE of interest is: $$u_t = -\frac{1}{2}(1-u)^2u_{x}-\delta^2u_{xx}+6\delta^2u(u_x)^2+3\delta^2u^2u_{xx}-u_{zz}+6u(u_z)^2+3u^2u_{zz}-\gamma u_{zzzz}$$ and the BCs are:

$$$$u_z=\lambda_0(1+u) \quad\text{on } z=0 \\ u_z=\lambda_H(1-u) \quad\text{on } z=1 \\ u_z-3u^2u_z+\gamma u_{zzz}=0 \quad\text{on } z=0 \\ u_z-3u^2u_z+\gamma u_{zzz}=0 \quad\text{on } z=1 \\$$$$


I want to apply the forward Euler Method on the linear terms in the PDE (handle them implicitly to reduce time step restrictions) and then handle the nonlinear terms explicitly after the fact. So after applying the discretized equation at the top, I would then solve for the nonlinear terms using: $$u_{j,k}^{n+1} = u_{j,k}^{n+1} + \Delta t N_{j,k}^n$$ where $$N_{j,k}^n$$ represents the discretized scheme for the nonlinear terms at time-level $$n$$ centered at the point $$(x_j,z_k)$$.

My real interest in asking this question is not so much to solve this particular equation with these specific BCs but to figure out a way to use some software to better handle this process of arranging a system of linear equations which stem from a single equation centered at different points into matrix-vector form so it can be solved using basic linear algebra techniques. Although in theory I know how to do this perfectly fine, it's tedious and susceptible to little, silly mistakes which are currently plaguing me.

• Why are you not simply using e.g. NDSolve to solve the PDE? Commented Feb 22, 2022 at 13:49
• I wanted to create my own Matlab script to solve this PDE. I have used NDSolve before and it's quite impressive but I'm hoping to develop my own methods that can deal with any subtleties as I add more intricacies and complications to the model Commented Feb 22, 2022 at 14:05
• @Mjoseph Could you show your PDE with boundary conditions? Commented Feb 22, 2022 at 14:40
• I added them as an EDIT to the original post - for some reason StackExchange wouldn't let me post it without editing the BCs as code. Commented Feb 22, 2022 at 15:49
• Check CoefficientArrays. It might do what you need. Commented Feb 22, 2022 at 16:28