I'm trying to code a matlab function that approximates a 1d diffusion equation with neumann boundary conditions, but I'm really struggling to understand the general approach for solving the problem. Could anyone tell me a general set of steps to approximating the equation using the theta method?
From what I understand it goes:
- create a matrix with size equivalent to the mesh size
- Fill the bottom row with values defined by the initial condition
- Fill the sides with values from the boundary conditions (how do i find these?)
- using the theta method formula below, calculate the mesh point above each initial value
$\frac{(u^{n+1}_j − u^n_j)}{∆t} = κ ((1 − θ)\frac{(u^n_j+1 − 2u^n_j + u^n_{j−1})}{h^2} + θ\frac{(u^n+1_j+1 − 2u^n+1_j + u^n+1_j−1)}{h^2})$
- repeat for each time step until the matrix is full
I hope thats at least close to the general process, but I cant understand how to work with the neumann boundary conditions.
I hope im not being too vague but any help would be great!!