Lets say we want to transmit a sequence of discrete data $\left\lbrace x[n] \right\rbrace$. But because we are living in analog world, the sequence must be modulated.
Call $T_s$ is symbol duration and use a set of orthonormal waveforms $\left\lbrace p_n(t) = p(t-nT_s), n \in \mathbb{Z} \right\rbrace$, (baseband) signal $x(t)$ can be written as
\begin{equation}
x(t) = \sum_{n} x[n] p(t - n T_s) \tag 1
\end{equation}
The equation (1) can be thought of modulating $x[n]$ to the n-th dimension of $x(t)$. To get back $x[n]$, we project $x(t)$ to this n-th dimension.
This projection is implemented by passing the received signal over the matched filter $q(t) = p^*(-t)$ and then sampling at $t=n T_s$. This projection is perfect if $g(t) = p(t) \star q(t)$ satisfies the the criterion $g(t) = 1$ if $t=0$ and $g(t) = 0$ if $t=m T_s\textrm{, } m \in \mathbb{Z}\textrm{, } m \neq 0$, i.e. Nyquist ISI criterion. Otherwise, the project will introduce error and by projecting to the n-th dimension, we don't get only $x[n]$ but also $x[m \neq n]$.
($\star$ means convolution)
This is bandlimited-caused ISI, not what you are asking but the model is necessary to understand the second ISI type, multipath-propagation-caused ISI.
Another source of ISI is multipath phenomenon. The received signal can be written as
\begin{equation}
y(t) = \sum_{i} a_i(t)x(t-\tau_i(t)) + w(t) \tag 2
\end{equation}
where $i$ is the index of physical path, $a_i(t)$ is its complex gain and $\tau_i(t)$ is the delay introduced by the path. $w(t)$ is AWG noise at receiver.
Apply the matched filter and sampling the output at $t = m T_s$,
\begin{eqnarray}
y[m] = y(t) \star q(t)|_{t = m T_s} = \sum_{i} a_i(m T_s) \sum_{n} x[n] g(m T_s - n T_s - \tau_i(m T_s)) + w[m] \\
y[m] \stackrel{l=m-n}{=} \sum_{l}x[m-l]h_l[m] + w[m]
\end{eqnarray}
with the channel tap $h_l[m]$, i.e. channel tap $l$ at the sample time $t=mT_s$,
\begin{equation}
h_l[m] = \sum_{i}a_i(m T_s)\times g(l T_s - \tau_i(m T_s)) \tag 3
\end{equation}
Now assume that the channel impulse reponse is invariant (a relaxed assumption could be that channel taps does not change during a coherence time, we call it underspread assumption),
$$y[m] = \sum_{l}x[m-l] \times h_l + w[m] \tag 4$$
The equation (4) is what you call FIR model. The equation (3) is the model of channel tap.
Remember Nyquist ISI criterion above ? $g(t) = 1$ if $t=0$ and $g(t) = 0$ if $t=m T_s\textrm{, } m \in \mathbb{Z}\textrm{, } m \neq 0$. Thus the power of the composite filter $g(t)$ must "concentrate" between $[-T_s/2, T_s/2]$ as in image for raised cosine pulse.
Thus the term $g(l T_s - \tau_i(m T_s))$ in Equation (3) is significant if and only if $-T_s/2 < |l \times T_s - \tau_i(m T_s)| < +T_s/2$.
We can interpret as the channel tap $l$ is contributed by the physical
path $i$ whose delay $\tau_i$ is around $t = l \times T_s$.
The channel impulse reponse can be thought being sampled by symbol duration $T_s$. And the number of tap depends on the relation between symbol duration $T_s$ and delay spread $\tau_m = \mathrm{max}_i (\tau_i)$, i.e. $L \approx \tau_m / T_s$.
As soon as $L > 1$, $y[m]$ is contributed by several $x[m], x[m-1], ...$ and it is multipath ISI. An equalizer will be needed to mitigate the effect of this ISI.
If $\forall \tau_i \ll T_s$, all physical path gains $a_i$ will contribute to only the channel tap $l=0$, we have flat fading, or single tap model.
What you are doing with matlab is playing with Equation (4) and you don't need to care about symbol duration and delay spread anymore.
xgenerated as shown in the code, where I don't know the sampling frequency. If I expand the above FIR model for the real domain case, say 3 channel coefficients, $h = [1,0.2,0.3]$ then for the $y(3) = h(1)*x(3) + h(2)*x(2) + h(3)*x(1)$ would be the output at time instant3. In theory it says that if the delay spread is greater that the symbol duration then we get ISI. Where is the ISI here, I don't understand. $\endgroup$