I'm afraid there's no shortcut. To get started you have to grab the bull by its horns. One has to read the manual, and get some simple example to work, and take it from there.

Please not that D02EBF has been phased out and is replaced by D02EJF.

The manual (short description) is online, at https://www.nag.co.uk/numeric/fl/manual/pdf/D02/d02ejf.pdf

It may take a while to gain experience in using these solvers, but the learning curve is fast, as soon as you have your first problem running. 

Thank you so much. I hope I can contact you in case of some doubts. 

Gustav Soederlind gave you already the pointer to the description of the routine. the bulk of work for you will be writing the subroutine f and pederv.  if you have already written simple f95 code - otherwise you anyway must first learn this at least on a basic level, and if your ode system isn't too complex, you also may first play with an equivalent routine which has almost all done this for you, such that only the body of fcn  must be done by yourself. I give you an simple example:

  y1 ' =  -3.0*y1 + y2 -x* y1*y2

 y2' = -y1- 10.0*y2 + x^2   

then the body of fcn would be

  f(1) = -3.0d0*y(1)+y(2)-x*y(1)*y(2)

 f(2) = -y(1)-10.0*y(2)+x**2 

in pederv you would code the partial derivatives of f,

jac(1,1)= -3.0d0 -x*y(2)

jac(1,2) = 1.d0 -x*y(1)

jac(2,1) = -1.d0

jac(2,2)= -10.0

small systems of this kind could be solved, writing only this body of fcn , using

http://numawww.mathematik.tu-darmstadt.de

you go in the section for ode's and choose there the implicit codes and there ''backward differentiation", "gears method".  you find there a description how this works and then in the "user defined" page the input forms for the body of fcn (no other work required) ,  the input form for the initial values of y and some choices for parameters, for example desired precision. hope this helps.