Let $$f_i(x)=\sum_{k=0}^{\infty}\frac{x^k}{(3k+i)!}$$
For $i\in \{0,1,2\}$.
First we manipulate the series for $f_1(x)$: $$f_1(x)=\sum_{k=0}^{\infty}\frac{x^k}{(3k+1)!}\implies f_1(x^3)=\sum_{k=0}^{\infty}\frac{x^{3k}}{(3k+1)!}\implies xf_1(x^3)=\sum_{k=0}^{\infty}\frac{x^{3k+1}}{(3k+1)!}$$
Taking the derivatives of both sides: $$\frac{d\left(xf_1(x^3)\right)}{dx}=\frac{d\left(\sum_{k=0}^{\infty}\frac{x^{3k+1}}{(3k+1)!}\right)}{dx}=\sum_{k=0}^{\infty}\frac{d\left(\frac{x^{3k+1}}{(3k+1)!}\right)}{dx}=\sum_{k=0}^{\infty}\frac{(3k+1)x^{3k}}{(3k+1)!}=\sum_{k=0}^{\infty}\frac{x^{3k}}{(3k)!}=f_0(x^3)$$
Performing similar manipulations on $f_2(x)$ yields $$\frac{d^2\left(x^2f_2(x^3)\right)}{dx^2}=f_0(x^3)$$
Lastly, if you consider the fact that each number is either congruent to either $0,1,$ or $2$ $\operatorname{mod} 3$, you get that $$e^x=\sum_{j=0}^{\infty}\frac{x^j}{j!}=\sum_{k=0}^{\infty}\frac{x^{3k}}{(3k)!}+\sum_{k=0}^{\infty}\frac{x^{3k+1}}{(3k+1)!}+\sum_{k=0}^{\infty}\frac{x^{3k+2}}{(3k+2)!}=f_0(x^3)+xf_1(x^3)+x^2f_2(x^3)$$
So $$e^x=f_0(x^3)+xf_1(x^3)+x^2f_2(x^3)$$
If we differentiate both sides wrt to $x$ twice we get $$\frac{d^2\left(e^x\right)}{dx^2}=\frac{d^2\left(f_0(x^3)\right)}{dx^2}+\frac{d^2\left(xf_1(x^3)\right)}{dx^2}+\frac{d^2\left(x^2f_2(x^3)\right)}{dx^2}$$
Using what we have already shown and the fact that $e^x$ is fixed under differentiation, we get $$e^x=\left(f_0(x^3)\right)''+\left(f_0(x^3)\right)'+f_0(x^3)$$
If we substitute $g(x)=f(x^3)$, we get $$e^x=g''(x)+g'(x)+g(x)$$ which is a simple ODE.
You can easily check that substituting $\frac{e^x}{3}$ for $g(x)$ will satisfy the above equation.
The general solution for the ODE is then given by $\frac{e^x}{3}+h(x)$ where $h(x)$ is the general solution for $$0=h''(x)+h'(x)+h(x)$$
From the basic theory of ODEs, we may find the genera $h(x)$ using the characteristic equation of the $h(x)$ ODE.
The roots are $$\frac{-1\pm i\sqrt{3}}{2}$$
Using this, the general solution for $g(x)$ becomes $\frac{e^x}{3}+Ae^{\frac{-1+ i\sqrt{3}}{2}x}+Be^{\frac{-1- i\sqrt{3}}{2}x}$ with $A$ and $B$ being arbitrary constants.
$f_0(x^3)$ is a special case $g(x)$ so by plugging in values of $x$ where $f_0(x^3)$ (or its first derivative) may be easily evaluated, we can solve for $A$ and $B$, and so solve for $f_0(x^3)$ and hence $f_0(x)$. Then using the relationships $f_1(x)$ and $f_2(x)$ have to $f_0(x)$, they may be solved for too.
Note: If there's anything wrong or confusing in my answer, feel free to edit it (and comment please) or just comment