Bruno,
In your code yP[0] == 0 is specified twice.
Linear systems of ODEs with constant coefficients are handled well by DSolve[], if your system is nonlinear consider using NDSolve (in case DSolve[] fails). Also, linear systems with constant coefficients can be solved explicitly. Here is an example:
(* X' = A.X + F *)
A = {{1, 0, 0, 0, 0}, {1/2, 1, 0, 0, 0}, {1/2, 1/2, 1, 0, 0}, {0, 0,
0, 1, 1/2}, {-1, 0, 0, 0, 1}} ;
X = {x1, x2, x3, x4, x5} ;
X0 = {1, 0, 0, 0, 0} ;
F = {t, Sin[t], 0, Cos[t], DiracDelta[t - 1]} ;
(* DSolve *)
Through[X[t]] /. DSolve[
{
Thread[D[Through[Flatten[X][t]], t] == A.Through[X[t]] + F],
Thread[Through[X[0]] == X0]
},
Through[X[t]],
t
] // Flatten // Simplify
(* Explicit formular *)
MatrixExp[t A].X0 +
MatrixExp[t A].Integrate[
Inverse[MatrixExp[x A]].(F /. t -> x), {x, 0, t},
Assumptions -> Element[t, Reals]] // Simplify
I.M.
