How to implement a dirac delta as input signal for this system
The output sequence y(n) is a linear regression of the input sequence x(n) which "keeps memory" of the last 2 samples of the input to compute the output.

Allocating the sequence x(n), required being a dirac pulse, into an array x[n] is a correct begin. Then you can develop the output sequence y(n) per input sequence x(n). For that I suggest you to exploit the functionality of matlab and to avoid the for loop:

y = zeros(1,20); % Allocate space for output
y(1) = b0*x(1); %initial condition
y(2:20) = b0*x(2:20) + b1*x(1:20); % for n &gt; 0


You should obtain an array of 20 samples of y(n) that you can easily plot using the function plot(x,y).
Note: since the dirac function is symbolic, you must decide by your own at which time the pulse should rise (e.g. at n=0 or otherwise). For that I suggest you the example in the documentation of the function called "Plot Dirac Delta Function".