Vikash linked derivations, so I'll capture the intuition. The Euler method numerically gets the first-derivative correct. RK4 gets the first four derivatives correct.

From the Taylor series, this means that RK4's first error term is from the 5th derivative term, which from the Taylor series has a (x-x_0)^5. = h^5 or if you prefer dx^5. So it's error is like the fifth power of h and is called 4th order. The Euler method is incorrect at the second derivative, which from the Taylor series has a (x-x_0)^2 = h^2. This is why it scales as h^2 and is called 1st order.

To understand what this really means, think about scaling your solution. You have a solution at h, and now you re-solve at h/2. How much did you change your error? For the Euler method your error is (h/2)^2 = h^2/4, so you have a quarter of the error from before. For the RK4 method you have (h/2)^5 = h^2/32 and so you have 1/32 the error from before! Thus you can see that as h gets smaller the higher order method gets better and better.

To better understand this, you should do the error analysis of the different methods which you can find in most of the standard books on numerical analysis. To summarize, if h is the step size, then local truncation error Euler's method is h^2 while for RK, 4th order it is h^5. The answer is essentially embedded in the formulation of the numerical schemes. There are even higher order RK methods which can provide even more accurate solutions.

http://www.math.unl.edu/~gledder1/Math447/EulerError

and

http://www.computer.org/csdl/proceedings/afips/1968/5071/00/50710467.pdf

also see the wikipedia article of both methods.

Great, I think I got it, thanks a lot 

Vikash linked derivations, so I'll capture the intuition. The Euler method numerically gets the first-derivative correct. RK4 gets the first four derivatives correct.

From the Taylor series, this means that RK4's first error term is from the 5th derivative term, which from the Taylor series has a (x-x_0)^5. = h^5 or if you prefer dx^5. So it's error is like the fifth power of h and is called 4th order. The Euler method is incorrect at the second derivative, which from the Taylor series has a (x-x_0)^2 = h^2. This is why it scales as h^2 and is called 1st order.

To understand what this really means, think about scaling your solution. You have a solution at h, and now you re-solve at h/2. How much did you change your error? For the Euler method your error is (h/2)^2 = h^2/4, so you have a quarter of the error from before. For the RK4 method you have (h/2)^5 = h^2/32 and so you have 1/32 the error from before! Thus you can see that as h gets smaller the higher order method gets better and better.

Dear Luciana

This book is useful in mathematica

www.springer.com/cda/content/.../cda.../9783540402459-c1.pdf?

http://www.math.bme.hu/~jtoth/Mma/Mathematica%20Cookbook.pdf

http://www.mathprogramming-intro.org/download/MathProgrammingIntro.pdf

Also, check out this textbook: http://comsizo.com.br/resolucoes/DennisZill_A%20First%20Course%20in%20Differential%20Equations%20with%20Modeling%20Applications.pdf

Thank you so much for your help. I am recently understanding how it works, and I need definitely those books and papers. 

No Problem

@ Luciana

The links I posted above should suffice. If there is still a problem, feel free to ask. Also, there are plethora of resources online!

The original question asks about accuracy, which was well discussed in the responses, but then goes on to raise the question of stability.   This is quite different and requires a separate analysis since it is possible to have a very accurate scheme which is unstable (especially for stiff problems) so the computed results deviate wildly from the true solution.

I would like to recommend the following book;

https://www.morebooks.de/store/gb/book/toward-absolutely-stable-numerical-scheme/isbn/978-3-659-68628-3

, which explains in detail a procedure to increase the numerical accuracy in the section 4.3 (4) in spite of a fully explicit scheme.  

Thank you so much for your contributions 

Dear, Thomas I. Seidman,

You are quite right! Simple example is numerical integration of function sqrt(x) from 0 to 1.

Luciana Chavez

The RK not only higher order, but it is the special polynomial approximation. There are many variants of RK of the same order. You also must know that if You solve the system the RK has some problem near the boundary.

Dear Drs. Thomas I. Seidman and Luciana Chavez

As you pointed out, numerical stability as well as accuracy is important. The RK schemes satisfy the possibilty condition for any time step interval Δt, which is a suffuicient condition to assure the numerical non-ossilation (unphysical non-oscillation) and boundedness condition due to accumulation of roundoff erros, which is explained in detail in the section 6.1.4 of Chaper 6 for one ODE and in the section 7.8.2 of Chapter 7 for advection-diffusion system eqations in my book as mentioned before:

https://www.morebooks.de/store/gb/book/toward-absolutely-stable-numerical-scheme/isbn/978-3-659-68628-3.

Dear Drs. Thomas I. Seidman and Luciana Chavez

I would like to correct one spelling "possibility" as "possitivity" in my tlast Answer as follows:

As you pointed out, numerical stability as well as accuracy is important. The RK schemes satisfy the possitivity condition for any time step interval Δt, which is a suffuicient condition to assure the numerical non-ossilation (unphysical non-oscillation) and boundedness condition due to accumulation of roundoff erros, which is explained in detail in the section 6.1.4 of Chaper 6 for one ODE and in the section 7.8.2 of Chapter 7 for advection-diffusion system eqations in my book as mentioned before: https://www.morebooks.de/store/gb/book/toward-absolutely-stable-numerical-scheme/isbn/978-3-659-68628-3.

the accuracy of any numerical method depends on (order of the method ) i.e. high order methods gives more accurate results

Thanks @ prof Sakai for the comprehensive answer.

There is a new method called piecewise analytic method (PAM).

Preprint Piecewise Analytic Method VS Runge-Kutta Method (Comparative Study)

Piecewise analytic method (PAM) is a new technique used for solving nonlinear differential equation. The advantages of PAM are( these points are proofed):

1. It gives a general analytic form that can be used in differentiation and integration.

2. It can solve highly non-linear differential equation.

3. The accuracy and error can be controlled according to our needs very easily.

4. It can solve problems which other famous techniques can’t solve.

5. In some cases, it gives the exact solution.