How to show a solution to an ODE system doesn't exist 
Show that each solution $(x(t), y(t))$ of the initial value problem
  $$
\left\{\begin{array}{cc}x' =&x^2+y \\ y' =&y^2+x \end{array}\right.\qquad
\left\{\begin{array}{cc}x(0) =&x_0 \\ y(0) =&y_0 \end{array}\right.
$$
  with $x_0>0$ and $y_0>0$ cannot exist on an interval of the form $[0,\infty)$.

I have been learning about different theorems to show existence but I am not sure how I would show this DNE. I put it into matlab using this code and got no solutions but is there a way to show this algebraically or some other way?
syms x(t) y(t) x0 y0

ode1 = diff(x) == x^2+y;
ode2 = diff(y) == y^2+x;
odes = [ode1; ode2]
S = dsolve(odes)

xSol(t) = S.x
ySol(t) = S.y

[xSol(t), ySol(t)] = dsolve(odes)

cond1 = x(0) == x0;
cond2 = y(0) == y0;
conds = [cond1; cond2];
[uSol(t), vSol(t)] = dsolve(odes,conds)

fplot(xSol)
hold on
fplot(ySol)
grid on
legend('xSol','ySol','Location','best')

 A: Firstly, the MATLAB program outputs the message

Warning: Explicit solution could not be found. 

Basically it means that an attempt to find a solution as a closed-form expression was unsuccessful. It doesn't mean that the solution does not exist; moreover,  for some value $\epsilon$ > 0, there exists a unique solution on the interval $t\in[-\epsilon,\epsilon]$ (by the Picard–Lindelöf theorem).
Consider the initial value problem
\begin{equation}\tag{1}
\left\{\begin{array}{lll}
\dot x&=&x^2\\
\dot y&=&y^2,\\
\end{array}
\right.\qquad
\left\{\begin{array}{l}
x(0)=x_0\\
y(0)=y_0,
\end{array}
\right.
\end{equation}
$x_0>0$, $y_0>0$. Its solution is
$$
x_{(1)}(t)=-\frac{1}{t-\frac{1}{x_0}},\qquad y_{(1)}(t)=-\frac{1}{t-\frac{1}{y_0}}.
$$
This solution does not exist on $[0,\infty)$ because it tends to $+\infty$ at $t\to1/x_0^{-}$ and $t\to1/y_0^{-}$.
The solution $x_{(2)}(t)$, $y_{(2)}(t)$ of the original initial value problem
\begin{equation}\tag{2}
\left\{\begin{array}{lll}
\dot x&=&x^2+y\\
\dot y&=&y^2+x,\\
\end{array}
\right.\qquad
\left\{\begin{array}{l}
x(0)=x_0\\
y(0)=y_0
\end{array}
\right.
\end{equation}
can not be expressed as a closed-form expression, but we know that in some interval $(0,T)$ (before it tends to +$\infty$)
$$
\dot x_{(2)}(t)=x^2+y>x^2=\dot x_{(1)}(t)>0,
$$
$$
\dot y_{(2)}(t)=y^2+x>y^2=\dot y_{(1)}(t)>0,
$$
** Update **
($x(t)>0$ and $y(t)>0$ because they are greater than the solution of the 
initial value problem $\dot x=y$, $\dot y= x$, $x(0)=x_0$, $y(0)=y_0$ which is
$x(t)= x_0\cosh t+y_0\sinh t$, $y(t)= x_0\sinh t+y_0\cosh t)$
** End of update **
thus,
$$
x_{(2)}(t)= x_0+\int_0^t \dot x_{(2)}(t)\,dt>x_0+\int_0^t \dot x_{(1)}(t)\,dt=
 x_{(1)}(t)=
-\frac{1}{t-\frac{1}{x_0}}
$$
$$
y_{(2)}(t)= y_0+\int_0^t \dot y_{(2)}(t)\,dt>y_0+\int_0^t \dot y_{(1)}(t)\,dt=
 y_{(1)}(t)=
-\frac{1}{t-\frac{1}{y_0}}.
$$
The components of the solution of the problem (2) are greater than the functions that tend to $+\infty$ in finite time, therefore, the solution doesn't exist on $[0,\infty)$.
A: Since the right hand sides of your system are even analytic functions of $x$ and $y$ it follows that for any initial point $(x_0,y_0)\in\bigl({\mathbb R}_{>0}\bigr)^2$there is an $h>0$ and a solution 
$$t\mapsto\bigl(x(t),y(t)\bigr)\quad (0\leq t<h),\qquad x(0)=x_0, \ y(0)=y_0$$
of the given IVP. But such a solution cannot live forever, because it explodes in finite time. 
Proof. By inspection one verifies that $x(t)>0$, $y(t)>0$ for all $t\geq0$, hence
$$x'(t)> x^2(t)\qquad(t\geq0)\ .$$
It follows that
$$\int_0^t{x'(\tau)\over x^2(\tau)}\>d\tau>\int_0^t 1\>d\tau=t\ ,$$
which expands to
$${1\over x_0}-{1\over x(t)}>t\ ,$$
or
$$x(t)>{1\over{1\over x_0}-t}\ .$$
This shows that the solution will explode before time $T:={1\over x_0}$.
A: Given the system 
$x' = x^2 + y,  \tag 1$
$y' = y^2 + x,  \tag 2$
we define the vector field
$\vec X(x, y) = \begin{bmatrix} x^2 + y \\ y^2 + x \end{bmatrix}, \tag 3$
which is defined on all of $\Bbb R^2$.  Consider the set
$\Omega = \{ (x, y) \in \Bbb R^2 \mid x \ge x_0, y \ge y_0 \} \subset \Bbb R^2; \tag 3$
the boundary $\partial \Omega$ of $\Omega$ is
$\partial \Omega = \{ (x, y_0) \mid x \ge x_0 \} \cup \{ (x_0, y) \mid y \ge y_0 \};  \tag 4$
we see that $\partial \Omega$ consists of two half-lines emanating from $(x_0, y_0)$, each one parallel to one of the $x$- and $y$-axes, and each extending indefinitely in the positive direction.
On $\partial \Omega$, we have
$x' = x^2 + y \ge x_0^2 + y_0 > 0, \tag 5$
and
$y' = y^2 + x \ge y_0^2 + x_0 > 0; \tag 6$
it follows from (5) and (6) that $X(x, y)$ points into $\Omega$ on $\partial \Omega$, and from this observation we conclude that any trajectory initialized on $\partial \Omega$ enters $\Omega$ immediately, and never leaves it;  that is, if $x(t_0), y(t_0) \in \partial \Omega$, then $(x(t), y(t)) \in \Omega$ for all $t > t_0$; likewise if $(x(t_0), y(t_0)) \in \Omega$, then $(x(t), y(t)) \in \Omega$ for all $ t > t_0$ as well.  
Now let $(x_1, y_1) \in \Omega$, and consider the integral curve $\gamma(t) = (x(t), y(t))$ where $(x(t_0), y(t_0)) = (x_1,  y_1)$; we have seen that $\gamma(t) \in \Omega$ for all $t \ge t_0$: thus,  on $\gamma(t)$, 
$x' = x^2 + y \ge x^2 + y_0, \tag 7$
from which
$\dfrac{x'}{x^2 + y_0} \ge 1; \tag 8$
thus, in a spirit similar to the arguments given by Christian Blatter and AVK, 
$\displaystyle \int_{t_0}^t \dfrac{x'(s)}{x^2(s) + y_0}ds \ge \int_{t_0}^t ds = (t - t_0); \tag 9$
we have
$\displaystyle \int_{t_0}^t \dfrac{x'(s)}{x^2(s) + y_0}ds = \dfrac{1}{\sqrt{y_0}}\tan^{-1} \dfrac{x(t)}{\sqrt{y_0}} - \dfrac{1}{\sqrt{y_0}}\tan^{-1} \dfrac{x(t_0)}{\sqrt{y_0}}$
$= \dfrac{1}{\sqrt{y_0}}(\tan^{-1} \dfrac{x(t)}{\sqrt{y_0}} - \tan^{-1} \dfrac{x(t_0)}{\sqrt{y_0}}), \tag{10}$
and so, via (9),
$\dfrac{1}{\sqrt{y_0}}(\tan^{-1} \dfrac{x(t)}{\sqrt{y_0}} - \tan^{-1} \dfrac{x(t_0)}{\sqrt{y_0}}) \ge t - t_0, \tag{11}$
and thus
$\tan^{-1} \dfrac{x(t)}{\sqrt{y_0}} - \tan^{-1} \dfrac{x(t_0)}{\sqrt{y_0}} \ge \sqrt{y_0}(t - t_0), \tag{12}$
whence
$\tan^{-1} \dfrac{x(t)}{\sqrt{y_0}} \ge  \sqrt{y_0}(t - t_0) + \tan^{-1} \dfrac{x(t_0)}{\sqrt{y_0}} \tag{13}$
and thus
$x(t) \ge \sqrt{y_0} \tan(\sqrt{y_0}(t - t_0) + \tan^{-1} \dfrac{x(t_0)}{\sqrt{y_0}}), \tag{14}$
the direction of the inequalities being preserved in (11)-(14) since $\tan$, $\tan^{-1}$, etc., are monotonically increasing functions of their respective arguments.
We note that for $x(t_0) \ge x_0$,
$0 < \tan^{-1} \dfrac{x(t_0)}{\sqrt{y_0}} < \dfrac{\pi}{2}, \tag{15}$
thus if we set
$t_{max} = t_0 + \dfrac{1}{\sqrt{y_0}}(\dfrac{\pi}{2} - \tan^{-1} \dfrac{x(t_0)}{\sqrt{y_0}}) \tag{16}$
we see that the solution $x(t)$ cannot possibly be continued as far as $t_{max}$, since from (14)
$x(t) \ge \sqrt{y_0} \tan(\sqrt{y_0}(t - t_0) + \tan^{-1} \dfrac{x(t_0)}{\sqrt{y_0}}) \to \infty \; \text{as} \; t \to t_{max}^-; \tag{17}$
the solution "blows up" in finite time.
