I have found an alternative and amazingly simple pattern for your sequences, that allows to calculate the generation function in a much more simple way.
Consider all blocks that end with $1$. There are two types of such blocks: $A$ which is either $1$ or $01$, and $B$ which is $000^*1$, i.e. at least two zeros followed by $1$ ($^*$ means repeat zero or more times).
Now, every sequence of $0$'s and $1's$ can be constructed as $\{A;B\}^*0^*$. That is we repeat blocks ending in $1$ and add any number of $0$'s at the end. The resulting sequence does not contain $X=00100$ iff two conditions are met: we don't have two consecutive blocks $B$, i.e. $BB$ is forbidden, and if we end the repetitive part with $B$ then the zero-tail cannot contain more than one $0$.
Overall,
$$A^*(BAA^*)^*\{B;B0;0^*\}$$
The generating functions (note, that if $Z$ has generating function $z(x)$, then $Z^*$ has generating function $\frac{1}{1-z(x)}$):
$$A=\{1;01\}:a(x)=x+x^2$$
$$B=000^*1:b(x)=\frac{x^3}{1-x}$$
And by the product lemma, for $A^*(BAA^*)^*\{B;B0;0^*\}$ we have
$$f(x)=\frac{1}{1-x-x^2}\frac{1}{1-\frac{x^4(1+x)}{(1-x)}\cdot\frac{1}{1-x-x^2}}\cdot\left(\frac{x^3}{1-x}+\frac{x^4}{1-x}+\frac{1}{1-x}\right)=$$
$$=\frac{1+x^3+x^4}{1-2x+x^3-x^4-x^5}.$$
Note, that $f(x)=1+2x+4x^2+8x^3+16x^4+31x^5+60x^6+116x^7+225x^8+437x^9+849x^{10}+\ldots$, and, indeed, we should have $n_l$ sequences of length $l$ where $n_l=2^l$ for $l\le 4$, $n_l=2^l-(l-4)2^{l-5}$ for $5\le l\le 7$, $n_l=2^l-(l-4)2^{l-5}+\sum_{k=1}^{l-7}k2^{k-1}$ for $8\le l\le 10$, etc.
HINT: You could start by finding a regular expression for the strings that contain the substring $bb$, irrespective of their length: the simplest is
$$(a+b)^*bb(a+b)^*\;.\tag{1}$$
Each of the subexpressions $(a+b)^*$ matches the set of all strings over $\{a,b\}$, so $(1)$ matches anything of the form $ubbv$ with $u,v\in\Sigma^*$, i.e., any string that has the substring $bb$ in it somewhere.
But we want only the strings $ubbv$ of odd length. Now $|uv|=|ubbv|-2$, so $ubbv$ has odd length if and only if $uv$ has odd length. And $|uv|=|u|+|v|$, so we want to generate all pairs of strings $u$ and $v$ such that $|u|+|v|$ is odd. That happens precisely when one of $|u|$ and $|v|$ is odd and the other is even.
Explain why the regular expression $$(aa+ab+ba+bb)^*\tag{2}$$ matches precisely the strings over $\Sigma$ that have even length.
Find a regular expression that matches precisely the strings over $\Sigma$ that have odd length; you may find it helpful to start with $(2)$ and modify or add to it.
Use the previous results to get a regular expression that matches precisely the words $ubbv\in\Sigma^*$ such that $|u|$ is odd and $|v|$ is even.
Do the same for words $ubbv$ such that $|u|$ is even and $|v|$ is odd.
Now put the pieces together to get a regular expression that matches precisely those words over $\Sigma$ that contain the substring $bb$ and have odd length.
Best Answer
For the particular example $M=2$ and $R=(10 | 1)^\ast 1^\ast$, I suppose if we first convert it to a CFG then there appears to be a general formula that involves convolutions. For example:
Then to get a string of length $n$, concatenate a string from $A$ of length $k_1$, a string from $C$ of length $k_2$, a string from $D$ of length $k_3$ and a string of length $A$ of length $k_4$, then you get
$N_E(n) = \sum_{k_1+k_2+k_3+k_4 = n} N_A(k_1)N_C(k_2)N_D(k_3) N_A(k_4)$
where $0 \leq k_i \leq n$. Ok, so there was no need to use a CFG, we could arrive here directly by identifying the parts!
Note $N_A(k_1) = 2^{k_1}$, $N_D(k_3) = 1$, and $N_E(k_4) = 2^{k_4}$. $N_C(k_2) = N_C(k_2-1) + N_C(k_2-2)$, Fibonacci recurrence with $N_C(1) = 1$ and $N_C(2) = 2$.
$N_{E,1}(n) = \sum_{k_1+k_2+k_3+k_4=n} 2^{k_1+k_4} F(k_3) = \sum_{a+b \leq n} 2^{a} F(b) = \sum_{b\leq n} (2^{n-b+1}-1) F(b)$. Looks like this can be attacked with generating functions.
This approach seems to generalize to other regular expressions. Oh and you said you know how to do it for exact matching. So for substring, just add an $A$ before and an $A$ after :) ?
Oh drats... but then you might be double-counting. You'd have to use inclusion-exclusion to subtract the number of strings where the pattern appears at least twice, add back where it appears at least three times, etc. But maybe it is okay.
For $ACDACDA$ (string appears at least twice): $N_{E,2} = \sum_{k_1+k_2+k_3+k_4+k_5+k_6+k_7=n} = N_A(k_1)N_C(k_2)N_A(k_4)N_C(k_5)N_A(k_7) = \sum_{a+b+c<n}(2^a-1) F(b)F(c)$.
Look for a general pattern?
$N_E(n) = N_{E,1}(n) - N_{E,2}(n) + N_{E,3}(n) + \ldots + N_{E,n}(n)$.
I didn't realize you had a generating function for strings matching $R$. Let us assume this now. Given the generating function $G_R(x)$ that contains number of strings of length $n$ matching $R$, let us figure out some way to obtain the generating function encoding the number of strings of length $n$ containing a substring that matches $R$.
So first, we do the same trick above and pad the sides with arbitrary strings. So the number of strings of length $n$ containing at least 1 substring matching $R$ is the sum $\sum_{a+b \leq n} M^a N_R(b) = \sum_{b \leq n} \frac{M^{n-b+1}-1}{M-1} N_R(b)$.
In general, the number of strings of length $n$ containing at least k substrings matching $R$ is the sum $\sum_{b_1+\ldots+b_k} \frac{M^{n-(b_1+\ldots+b_k)+1}-1}{M-1} N_R(b_1)\cdots N_R(b_k)$.
As a generating function, let $H_k(x)$ correspond to containing at least $k$ substrings.
$H_k(x) = \sum_{j} N_{R,k}(j) x^j = \sum_{j} \sum_{a+b_1+\ldots+b_k = j} C(j-a) x^{j-a} N_R(b_1) x^{b_1} \cdots N_R(b_k) x^{b_k}$
This is a $k$-fold convolution between the generating function of all $M$-strings, and $k-1$ copies of $G_R$.
Now we encode inclusion exclusion:
$H(x) = \sum_{j} \sum_{i=1}^j (-1)^{i-1} N_{R,\ i}(j) x^j = \sum_i (-1)^{i-1} \sum_{j \geq i} N_{R,\ i}(j) x^j $, and this is just the alternating sum of the corresponding partial generating functions. So... I'm going to stop here for the moment.