Finite Difference SolutionsOne way to solve partial differential equations is to form a spatial discretization and solve them with finite differences. This example considers a basic finite difference solution for the Poisson equation on a unit square. This sets up the problem for an 8080 grid with . An abbreviated form of the variable matrix is printed out. In[1]:= 
Out[3]//Short=

This is the finite difference formula for the differential operator. In[4]:= 
In this step the application of the finite difference formula to the variable matrix is made. The boundary conditions are incorporated by allowing an overhang of one variable with a coefficient of zero around the boundary. In[5]:= 
Here equations for the solutions on the grid are formed. In[6]:= 
Out[6]//Short=

The approximate solution to the Poisson equation will be found by solving these equations. One way to solve the equations is to use Solve, which will recognize the equations as linear and sparse and accordingly use a sparse matrix solver. Another method that also gives you the possibility of solving additional Poisson problems is to construct a matrix from the equations. The matrix represents the linear operator when the grid is flattened into a vector. In[7]:= 
Out[7]=

Now the solution is quickly found using LinearSolve. In[8]:= 
Out[8]=

This repartitions the solution to restore the original structure. In[9]:= 
This makes a contour plot of the solution. In[10]:= 
If you wanted to solve several Poisson equations with the same geometries but different righthand functions, you could compute a factorization of the matrix and use it repeatedly. The way that Mathematica allows you to repeat the factorization is described in the section Saving the Factorization. In[11]:= 
Out[11]=

Now you just need to give it a vector representing the function on the grid. In[12]:= 
This passes the vector to the factorization, thereby finding the solution. In[14]:= 
It should be noted that this solution is only available at discrete points in the interior of the region. It does not include the value on the boundary or at other points in the interior. If the solution was to be used for other computational purposes this might be a limitation. PadLeft and PadRight could be used to add the boundary value and Interpolation used to generate an interpolating function for the region. In[15]:= 
Out[17]=

This uses Plot3D to plot the solution. In[18]:= 
Now the solution can be computed at any position in the region. In[19]:= 
Out[19]=

