ODE solver giving unexpected result I have the following coupled ODEs:
$$
\begin{cases}
F'(s) & = \lambda F(s) - \lambda (1 - (1 - F(s))^2) + \lambda H'(s)\\
H'(s) & = 1 - H(s) - (1 - F(s))^2,
\end{cases}
$$
with boundary condition $F(0) = 1-\lambda, H(0) = 0$, here $\lambda \in ]0,1[$. The solution for $F$ I am looking for is a Cummulative Distribution Function (CDF), but when I solve these coupled ODEs with these boundary conditions, I get a result which is not a CDF at all. The question is: how to find the correct solution OR did I make a mistake in deriving these ODEs and do these ODEs indeed not give a CDF as a possible solution?
I am using the Python solve odeint
 A: You can rewrite the equation in the following simple formulation (by simply inserting the definition of $H'$):
\begin{align}
F' & = \lambda (F - H) \\
H' & = 1 - H - (1 - F)^2
\end{align}
This equation hast two fixed points:
\begin{align}
H &= 0 , \ F = 0 \\
H &= 1, \ F = 1
\end{align}
So you can imagine that if you start somewhere between $0$ and $1$ you might converge to one of these fixed points. But is one attractive and one repulsive?
To answer this question we actually look at another property. Since $F$ is a CDF it's supposed to be monotonically increasing (and converging to $1$). So we need $$D =  H - F < 0$$
So we write the differential equation for $D$ and $F$:
\begin{align}
F' &= -\lambda D \\
D' &= (\lambda - 1)D + F(1-F)
\end{align}
We see that there is no particular reason for the latter equation to stay negative. In any case we can solve the equation numerically. First we see that the negativity of $D$ is not satisfied. Furthermore on large scales we converge to $0$ and not to $1$.

A: In the integral equation as given in the comments
$$
F(s) = F(0) + \lambda \int_0^s F(u) du + \lambda \int_0^s(1-(1-F(u))^k) (1-e^{u-s}) du
$$
set again $H(s)= \int_0^s(1-(1-F(u))^k)e^{u-s}du$. Then
$$
H'(s)+H(s)=(1-(1-F(s))^k)
$$
Now taking the derivative of the first equation gives
$$
F'(s)=λ(F(s)+(1-(1-F(u))^k)-H'(s))=λ(F(s)+H(s)).
$$
The stationary points of this system satisfy $0=f+h$, $h=1-(1-f)^k$ which for $k=2$ results in  the quadratic equation $$0=1+f-(1-f)^2=3f-f^2$$ which has a stationary point at $f=0$ and $f=3$. Numerical integration confirms that result.
Now if the sign of the first integral were changed, 
$$
F(s) = F(0) - \lambda \int_0^s F(u) du + \lambda \int_0^s(1-(1-F(u))^k) (1-e^{u-s}) du
$$
then the system 
\begin{align}
F(s)&=λ(-F(s)+H(s))\\
H'(s)&=1-H(s)-(1-F(s))^k
\end{align}
has its stationary points at $f=h=1-(1-f)^k$, that is, at $f=h=0$ and $f=h=1$, as required for a CDF that goes from asymptotically zero to asymptotically one. One would have to explore the stability and basins of attraction for these stationary points, monotonicity is not a given as in the case of scalar first order ODE.
A numerical example integration gives 
def deriv(y,a): F,H = y; return [a*(-F+H), 1-H-(1-F)**2]
y0 = [0.3, 0]
a = 0.2
tp = np.linspace(0,30,151)
tm = np.linspace(0,-5,21)
solp=odeint(lambda y,t: deriv(y,0.2), [ 0.1, 0], tp)
solm=odeint(lambda y,t: deriv(y,0.2), [ 0.1, 0], tm)
plt.plot(tp,solp[:,0],tm,solm[:,0]); plt.grid(); plt.show()

with the result

However, at $s=0$ with $H(0)=0$ one has $F'(0)=-λF(0)$ so that the solution has an initial negative slope which is untypical for CDF.
