@Mariusz Iwaniuk  thanks for your help. but i still have some issues with the code i develop.

                              NULL
> restart;
> dsys := {diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-H[a]*x[3](y)-V[0]*x[3](y)-B[r]*x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -H[a]*P[m]*x[3](y)-V[0]*P[m]*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -N[t]*x[7](y)^2-V[0]*P[r]*x[7](y)-N[b]*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -V[0]*S[c]*x[9](y)+N[t]^2*x[7](y)^2/N[b]+N[t]*P[r]*x[7](y)/N[b]+N[t]*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 0, x[4](0) = 0, x[5](0) = 0, x[6](0) = 0, x[7](0) = -1, x[8](0) = 1, x[9](0) = -.68};
 /                                                              
 | d                d                      d                    
< --- x[1](y) = 1, --- x[2](y) = x[3](y), --- x[3](y) = -x[6](y)
 | dy               dy                     dy                   
 \                                                              

   - H[a] x[3](y) - V[0] x[3](y) - B[r] x[8](y), 

   d                     
  --- x[4](y) = x[5](y), 
   dy                    

   d                                                    
  --- x[5](y) = -H[a] P[m] x[3](y) - V[0] P[m] x[5](y), 
   dy                                                   

   d                      d                         2
  --- x[6](y) = x[7](y), --- x[7](y) = -N[t] x[7](y) 
   dy                     dy                         

   - V[0] P[r] x[7](y) - N[b] x[7](y) x[9](y), 

   d                      d                              
  --- x[8](y) = x[9](y), --- x[9](y) = -V[0] S[c] x[9](y)
   dy                     dy                             

         2        2                                             
     N[t]  x[7](y)    N[t] P[r] x[7](y)                         
   + -------------- + ----------------- + N[t] x[7](y) x[9](y), 
          N[b]              N[b]                                

  x[1](0) = 0, x[2](0) = 0, x[3](0) = 0, x[4](0) = 0, 

  x[5](0) = 0, x[6](0) = 0, x[7](0) = -1, x[8](0) = 1, 

                 \ 
                 | 
  x[9](0) = -0.68 >
                 | 
                 / 
> sys1 := eval(dsys, {B[r] = 1, H[a] = 5, N[b] = .1, N[t] = .1, P[m] = .8, P[r] = 10, S[c] = 1, V[0] = 1});
 / d                d                     
{ --- x[1](y) = 1, --- x[2](y) = x[3](y), 
 \ dy               dy                    

   d                                            
  --- x[3](y) = -x[6](y) - 6 x[3](y) - x[8](y), 
   dy                                           

   d                     
  --- x[4](y) = x[5](y), 
   dy                    

   d                                        
  --- x[5](y) = -4.0 x[3](y) - 0.8 x[5](y), 
   dy                                       

   d                     
  --- x[6](y) = x[7](y), 
   dy                    

   d                        2                                     
  --- x[7](y) = -0.1 x[7](y)  - 10 x[7](y) - 0.1 x[7](y) x[9](y), 
   dy                                                             

   d                      d                    
  --- x[8](y) = x[9](y), --- x[9](y) = -x[9](y)
   dy                     dy                   

                         2                      
   + 0.1000000000 x[7](y)  + 10.00000000 x[7](y)

   + 0.1 x[7](y) x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 0, 

  x[4](0) = 0, x[5](0) = 0, x[6](0) = 0, x[7](0) = -1, 

                              \ 
  x[8](0) = 1, x[9](0) = -0.68 }
                              / 
> sol := dsolve(dsys, numeric, method = classical[rk4], stepsize = 0.1e-1);
Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)
proc(x_classical)  ...  end;
> plots:-odeplot(sol, [[y, x[8](y)]], y = 0 .. 2, color = ["Red"]);
Warning, could not obtain numerical solution at all points, plot may be incomplete