The substitution method is probably best for what you are trying to do.
On the other hand, if you want to compute the derivatives of $e^{-x^2}$, then you will have to use a combination of the chain rule and the product rule.
$$f(x) = e^{-x^2}$$
$$f'(x) = -2x e^{-x^2}$$
$$f''(x) = - 2 e^{x^2}+ (-2x)^2 e^{-x^2} = (4x^2 -2) e^{-x^2}$$
$$f'''(x) = 8x e^{-x^2} - (4x^2 - 2) (2x) e^{-x^2} = (-8x^3 + 12x) e^{-x^2}$$
This yields:
$f(0)=1$, $f'(0) = 0$, $f''(0) = -2$, $f'''(0) = 0$, etc.
This will be rather difficult.
Let's first of all suppose there are no problems, i.e. the functions you are interested in are infinitely differentiable and continuous in the range of values you are interested in.
Then, even with the simple example that you know the Maclaurin series of $f(x)$ and want the Taylor series of $f(x)$ about another point $a$, this is a cumbersome endeavour. Here is the procedure.
With the known Maclaurin series (you memorized the $m_k$)
$f(x) = \sum\limits_{k = 0}^{\infty} \frac{f^{(k)} (0)}{k!} x^k = \sum\limits_{k = 0}^{\infty} m_k x^k$
you require the Taylor series
$f_a(x) = \sum\limits_{k = 0}^{\infty} c_k (x-a)^k$
Now the coefficients are given by
$c_k = \frac{f^{(k)} (a)}{k!}$
which in turn you can get by differentiating your MacLaurin series and evaluate at $a$:
$c_k = \frac{f^{(k)} (a)}{k!} = \sum\limits_{n = k}^{\infty} \frac{n!}{k! (n-k)!} m_n a^{n-k}$
E.g.
$c_0 = \sum\limits_{n = 0}^{\infty} m_n a^{n}$
So you need to sum an infinite series to arrive at your new coefficients - big labour.
Another method is to identify the coefficients of all powers of $(x-a)$. To do so, write the MacLaurin series
$f(x) = \sum\limits_{k = 0}^{\infty} m_k ((x-a)+a)^k
= \sum\limits_{k = 0}^{\infty} m_k \sum_{j=0}^k {k \choose j} (x-a)^j (a)^{k-j}$
This leads to the same result.
Things get worse if you consider $g(f(x))$.
I guess this already answers your question: it's not a smooth technique to calculate the Taylor expansion from the known MacLaurin series.
Best Answer
We can form a series that is somewhat easier to work with by noting that $ \ \frac{d}{dx} \ [e^{x^3}] \ = \ 3x^2 · e^{x^3} \ \ , $ so that we produce the anti-derivative series $$ G(x) \ \ = \ \ \frac13 \ e^{x^3} \ \ = \ \frac13 \ \sum_{n \ = \ 0}^{\infty} \ \frac{1}{n!} \ (x^3)^n \ \ = \ \ \frac13 \ \sum_{n \ = \ 0}^{\infty} \ \frac{x^{3n}}{n!} \ \ , $$
with $ \ G'(x) \ = \ g(x) \ = \ x^2 e^{x^3} \ . $ The series for $ \ G(x) \ $ is non-zero only for every third term, that is to say,
$$ G^{(n)}(0) \ \ = \ \ \left\{ \begin{array}{rcl} a_{n} \ \ , & \ \text{for} \ n = 3k \ \ , \ \ k \ \ \text{a non-negative integer} \\ 0 \ \ , & \text{otherwise} \end{array} \right. \ \ . $$
Since $ \ g^{(40)}(0) \ = \ G^{(41)}(0) \ \ , $ and $ \ 41 \ $ is not a multiple of $ \ 3 \ \ , $ we can simply conclude that $ \ g^{(40)}(0) \ = \ 0 \ \ . $ (Had the problem asked for $ \ g^{(41)}(0) \ \ \text{or} \ \ g^{(50)}(0) \ \ $ (as appears in one of the answers -- the original version of the post?), we would have more work to do, as discussed below.
$$ \ \ $$
Another way we can approach the problem is to make use of the "higher-derivatives" version of the Product Rule, $$ \frac{d^n}{dx^n} \ [f_1(x) · f_2(x)] \ \ = \ \ \sum_{k \ = \ 0}^{n} \ \binom{n}{k} \ f_1^{(k)}(x) · f_2^{(n-k)}(x) \ \ , $$
which behaves like a binomial expansion. This will be of help to us if we define, as we did above, $ \ G(x) \ = \ \frac{1}{3} · e^{x^3} \ = \ \frac{1}{3} · e^u \ = \ y \ \ , $ so that $ \ g(x) \ = \ G'(x) \ = \ u' · y \ \ . $ We will then set $ \ f_1 = u' \ $ and $ \ f_2 = y \ \ $ in making our calculations.
This may look daunting for the derivative we wish to determine, but there are features of the function that lead to considerable simplification. We have $ \ G(0) = \frac13 \ $ and $ \ u = x^3 \ \Rightarrow \ u' = 3x^2 \ , \ u'' = 6x \ , \ u''' = 6 \ \ , $ and $ \ u^{(k)} = 0 \ $ for $ \ k \ge 4 \ \ . $ What this tells us is that the number of terms in the expansion for successive derivatives is quite limited. Moreover, since we are evaluating the derivatives of $ \ g(x) \ $ at $ \ x = 0 \ $ , the only non-zero value for the $ \ u^{(n)} $ 's is $ \ u''' = 6 \ \ . $
With this in mind, we see that
$ g(0) \ = \ G'(0) \ = \ y'(0) \ = \ 0 · \frac13 \ = \ 0 \ \ , $
$ g'(x) \ = \ G''(x) \ = \ y'' \ = \ u''· y \ + \ u'· y' \ \ \Rightarrow \ \ g'(0) \ = \ 0 · \frac13 \ + \ 0 · 0 \ = \ 0 \ \ , $
$ g''(x) \ = \ y''' \ = \ u'''· y \ + \ 2 · u''· y' \ + \ u'·y'' \ \ \Rightarrow \ \ g''(0) \ = \ 6 · \frac13 \ + \ 2 · 0 · 0 \ + \ 0·0 \ = \ 2 \ \ . $
Here, we attain our first non-zero derivative with $ \ G'''(0) \ = \ g''(0) \ \ . $ With the next derivative, our expressions begin to include derivatives of $ \ u \ $ which are simply zero.
$ g'''(x) \ = \ y^{(4)} \ = \ u^{(4)}· y \ + \ 3 · u'''· y' \ + \ 3·u''·y'' \ + \ u'·y''' $ $ \Rightarrow \ \ g'''(0) \ = \ 0 · \frac13 \ + \ 3 · \frac13 · 0 \ + \ 3·0·0 \ + \ 0·6 \ = \ 0 \ \ , $
$ g^{(4)}(x) \ = \ y^{(5)} \ = \ u^{(5)}· y \ + \ 4 · u^{(4)}· y' \ + \ 6·u'''·y'' \ + \ 4·u''·y''' \ + \ u'·y^{(4)} $ $ \Rightarrow \ \ g^{(4)}(0) \ = \ 0 \ + \ 0 \ + \ 6 · \frac13 · 0 \ + \ 4·0·2 \ + \ 0 \ = \ 0 \ \ , $
and
$ g^{(5)}(x) \ = \ y^{(6)} \ = \ u^{(6)}· y \ + \ 5 · u^{(5)}· y' \ + \ 10·u^{(4)}·y'' \ + \ 10·u'''·y''' \ + \ 5·u''·y^{(4)} \ + \ u'·y^{(5)} $ $ \Rightarrow \ \ g^{(5)}(0) \ = \ 0 \ + \ 0 \ + \ 0 \ + \ 10·6·2 \ + \ 0 \ + \ 0 \ = \ \binom{5}{3}·6·2 \ = \ 120 \ \ . $
What becomes clear from this examples is that the non-zero derivatives for non-negative integer $ \ k \ $ are given by the recursion relation
$$ g^{(3k + 2)}(0) \ = \ y^{(3k + 3)}(0) \ = \ \binom{3k+2}{3k}·u'''·y^{(3k)}(0) \ = \ \binom{3k+2}{3k}·6·y^{(3k)}(0) \ \ , \ \ y(0) \ = \ \frac13 \ \ . $$
Hence,
$$ g^{(8)}(0) \ = \ y^{(9)}(0) \ = \ \binom{8}{6}·6·y^{(6)}(0) \ = \ \binom{8}{6}·6· \left[\binom{5}{3}·6·2 \right] \ = \ 20,160 \ \ , $$ $$ g^{(11)}(0) \ = \ y^{(12)}(0) \ = \ \binom{11}{9}·6·y^{(9)}(0) \ = \ 55·6·20,160 \ = \ 6,652,800 \ \ , $$
and so forth. We observe that while $ \ g^{(41)}(0) \ = \ y^{(42)}(0) \ $ would be enormous, $ \ g^{(40)}(0) \ = \ y^{(41)}(0) \ $ is not in this set of derivatives, and so is equal to zero.
$$ \ \ $$
We might like to have something a bit simpler to work with than the recursion relation for the non-zero derivatives. If we write out the last example in full, we obtain
$$ g^{(11)}(0) \ \ = \ \ \binom{11}{9}·6·\binom{8}{6}·6· \binom{5}{3}·6·\binom{2}{0}·6·\frac13 $$ $$ = \ \ \frac{11·10·9·8·7·6·5·4·3·2·1}{2 \ · \ 9 \ · \ 2 \ · \ 6 \ · \ 2 \ · \ 3 \ · \ 2} \ · \ 6^3 \ · \ 2 \ \ = \ \ \frac{11!}{2^3 \ · \ 3^3 \ · \ (3·2·1) \ · \ 2} \ · \ 6^3 \ · \ 2 \ \ = \ \ \frac{11!}{3!} \ \ . $$
The general calculation yields $ g^{(3k + 2)}(0) \ = \ \frac{(3k+2)!}{k!} \ \ , $ for integer $ \ k \ge 0 \ $ ; this completes the derivation characterization given above
$$ G^{(n+3)}(0) \ \ = \ \ \left\{ \begin{array}{rcl} \frac{(3k+2)!}{k!} \ \ , & \ \text{for} \ n = 3k \ \ , \ \ k \ \ \text{a non-negative integer} \\ 0 \ \ , & \text{otherwise} \end{array} \right. \ \ . $$
Inserting this into the general term for a Maclaurin series, we then produce
$$ g(x) \ \ = \ \ x^2 · e^{x^3} \ \ = \ \ \sum_{m \ = \ 0}^{\infty} \ \frac{g^{(m)}(0) }{m!} \ x^m $$ [removing the "zero terms"] $$ = \ \ \sum_{k \ = \ 0}^{\infty} \ \frac{(3k+2)! \ / \ k! }{(3k+2)!} \ x^{3k+2} \ \ = \ \ \sum_{k \ = \ 0}^{\infty} \ \frac{x^{3k+2} }{k!} \ \ = \ \ x^2 \ + \ x^5 \ + \ \frac12 · x^8 \ + \ \frac16 · x^{11} \ + \ \ldots \ \ , $$
in agreement with the method of substitution into the series for $ \ e^u \ $ and Gary's comment.