filter() function implements a difference equation by calculating
for every sample index $n$, where $y[n]$ is the output sequence, and $x[n]$ is the input sequence.
The coefficients $a_i$ and $b_i$ are stored in the vectors
b. Note that normally the first element of the vector
a is chosen as $a_0=1$. From (1) you see that in order to implement the difference equation you need to store past input and output values. This is usually done with a circular buffer.
Normally you don't talk about "designing a difference equation", but you design a filter, i.e. you compute the coefficients $a_i$ and $b_i$ in such a way that the filter performs some chosen task. There is a large variety of filter design algorithms available. Take a look at this answer for some basic background information on filter design (this is a huge topic in itself!).
In most cases you would use some filter design package (e.g. available in Matlab or Octave), design the filter (i.e. compute the filter coefficients), and then you can implement the filter according to Eq. (1). Writing your own filter design routine can be a nice exercise, but since filter design is a very mature topic, you'll basically be re-inventing the wheel. Implementing the filter is a different thing, and it can be useful to write your own implementation if you want to perform some real-time filtering.
Note that a direct implementation of Eq. (1) is only one of many ways to implement a filter. There are many possible implementations, all with their advantages and disadvantages. They all implement Eq. (1), but they perform computations in different orders, and they use different intermediate variables, which results in very different numerical behaviors. Of course, they would all be equivalent if it were possible to perform exact computations without any numerical errors. Especially for fixed-point implementations it is very important to choose an appropriate structure, because a direct implementation of Eq. (1) will be disastrous for larger filter orders (i.e. larger values of $M$ and $N$).
Note that Eq. (1) is a time-domain implementation. If you're interested in frequency domain implementations, search for "overlap-add" and "overlap-save" methods. However, I would start with implementing a time-domain filtering routine if you're working on your first implementation. Also note that for frequency domain filtering you need an FFT routine.
If you can use floating point arithmetic, I would as a first step try to implement Eq. (1) directly, and compare your results with Matlab's
filter() function. For computing the coefficients you should use Matlab's filter design functions. If you then have any more specific questions you could come back and formulate a new question.