Let $M$ be a smooth Manifold.
Let $p \in M$.
Let $U$ be an open neighborhood of $p$.
We denote by $C^{\infty}(U)$ the set of real valued smooth functions on $U$.
Let $\Lambda_p = \bigcup C^{\infty}(U)$, where $U$ runs through all open neighborhoods of $p$.
Let $f, g \in \Lambda_p$.
Suppose $f \in C^{\infty}(U)$ and $g \in C^{\infty}(V)$.
If there exists an open neighborhood $W$ of $p$ such that $W \subset U \cap V$ and $f|W = g|W$,
we say $f$ and $g$ are equivalent.
This is an equivalence relation on $\Lambda_p$.
We denote by $\mathcal{O}_p$ the set of equivalence classes on $\Lambda_p$.
Clearly $\mathcal{O}_p$ is an $\mathbb{R}$-algebra.
Let $f \in C^{\infty}(U)$, where $U$ is an open neighborhodd of $p$. We denote by $[f]$ the equivalence class containing $f$.
A derivation of $\mathcal{O}_p$ is a linear map $D\colon \mathcal{O}_p \rightarrow \mathbb{R}$ such that
$D(fg) = D(f)g(p) + f(p)D(g)$ for $f, g \in \mathcal{O}_p$.
The set $T_p(M)$ of derivations of $\mathcal{O}_p$ is a vector space over $\mathbb{R}$ and is called the tangent space at $p$.
Let $\epsilon > 0$ be a positive real number.
We denote by $\Gamma_p(\epsilon)$ the set of smooth curves $\gamma \colon (-\epsilon,\epsilon) \rightarrow M$ such that $\gamma(0) = p$.
Let $\Gamma_p = \bigcup_{\epsilon>0} \Gamma_p(\epsilon)$.
Let $(U, \phi)$ be a chart such that $p \in U$.
Let $\gamma_1, \gamma_2 \in \Gamma_p$.
Then $\gamma_1$ and $\gamma_2$ are called equivalent at $0$ if $(\phi\circ\gamma_1)'(0) = (\phi\circ\gamma_2)'(0)$. This definition does not depend on the choice of the chart $(U, \phi)$.
This defines an equivalence relation on $\Gamma_p(M)$.
Let $S_p(M)$ be the set of equivalence classes on $\Gamma_p(M)$.
For $\gamma \in \Gamma_p(M)$, we denote by $[\gamma]$ the equivalence class containing $\gamma$.
We will define a map $\Phi\colon S_p(M) \rightarrow T_p(M)$.
Let $c \in S_p(M)$.
Choose $\gamma \in \Gamma_p(M)$ such that $c = [\gamma]$.
Let $f \in C^{\infty}(U)$, where $U$ is an open neighborhood of $p$.
We write $D_c([f]) = (f\circ\gamma)'(0)$ for $f \in C^{\infty}(U)$.
Clearly $D_c$ is well defined and does not depend on the choice of $\gamma$.
Clearly $D_c \in T_p(M)$.
Hence we get a map $\Phi\colon S_p(M) \rightarrow T_p(M)$ such that $\Phi(c) = D_c$.
We claim that $\Phi$ is bijective.
Let $c, e \in S_p(M)$.
Suppose $D_c = D_e$.
Suppose $c = [\gamma]$ and $e = [\lambda]$.
Let $(U, \phi)$ be a chart such that $p \in U$.
Let $\pi_i:\mathbb{R}^n \rightarrow \mathbb{R}$ be the $i$-th projection map: $\pi_i(x_1,\dots,x_n) = x_i$.
We denote by $\phi^i$ by $\pi_i\circ\phi$.
Since $D_c([\phi^i]) = D_e([\phi^i])$, $(\phi^i\circ\gamma)'(0) = (\phi^i\circ\lambda)'(0)$.
Hence $(\phi\circ\gamma)'(0) = (\phi\circ\lambda)'(0)$.
Hence $\gamma$ and $\lambda$ is equivalent.
Thus $\Phi$ is injective.
Let $D \in T_p(M)$.
Let $(U, \phi)$ be a chart such that $p \in U$.
We assume that $\phi(p) = 0$.
We define $\phi^i$ for $i = 1,\dots,n$ as above.
Let $D([\phi^i]) = a_i$ for $i = 1,\dots,n$.
There exists $\epsilon > 0$ such that $(a_1t,\dots,a_nt) \in \phi(U)$ for every $t \in (-\epsilon, \epsilon)$.
Let $\gamma(t) = \phi^{-1}(a_1t,\dots,a_nt)$ for $t \in (-\epsilon, \epsilon)$.
Then it's easy to see that $\Phi([\gamma]) = D$.
Hence $\Phi$ is surjective and we are done.
The prototypical example for this situation is something like the vector field $-x \partial_x$ on a one-dimensional manifold. This has a zero at $x=0$, and indeed the flow at $x=0$ is constant; but this does not stop the flow starting at any other point from being maximally defined.
Your mistake is in confusing "defined on all of $\mathbb R$" with a statement about the image of the curve. In my example above, the integral curve starting at $x_0$ is simply $x(t) = x_0 e^{-t}$, which is defined on all of $\mathbb R$, but never hits $0$ - it just asymptotically approaches it. The fact that the image of the curve could be extended (i.e. the curve could be extended after a reparametrization) does not mean that the curve itself is not maximally defined.
Best Answer
This capability is implemented in many programs, including matlab, mathematica, and matplotlib.
If you want to know the underlying algorithms, they work by numerically solving $$\frac{d}{dt} \mathbf{x}(t)=\mathbf{v}(\mathbf{x}(t))$$ where $\mathbf{v}$ is your vector field.