Improved Euler method for second order ODE I am trying to solve the simple harmonic oscillator problem with various Euler methods. Having managed to solve it with simple and modified Euler methods now I am trying to solve it with the improved Euler method which is a bit like Runge-Kutta. My concern is that when I plot the graph of energy vs. time it oscillates though it also increases at the same time. Since x(position) is actually the integral of v(velocity) is it possible to solve it numerically?
The equations I have used are as follows:
$$ \frac{dv}{dt}= a= \frac{F}{m}=-\frac{kx}{m}$$
and the following first-order equation for the position:
$$ \frac{dx}{dt}= v $$
EDIT: The initial conditions and the constants are
$$x(0) = 0, v(0) = 5,  \ \ k = m = 1,\ h=0.1, \ N=201  $$
Yet my concern is that these equations are not independent and I think the improved Euler method cannot be applied. Do you have any suggestions on how to approach the problem? I have also added the resulting graph to illustrate my point. I have updated position and velocity as follows and in the depicted order:
$$ x_{0}=x \\ x = x + hv \\ v = v - h(\frac{kx}{m} + \frac{kx_{0}}{m})/2 \\ E = (mv^2+kx^2)/2$$

There must be a mistake in this implementation and I would be grateful if someone point it out or makes a suggestion.
Thanks in advance
EDIT: Latest C Implementation
for(i = 0; i < N ; i++)
    {
        fprintf(fp, "%f \t %f \t %f \t %f \n", x, v, E, t);
        xint = x + (h) * v;
        vint = v - (h) * k * x / m;
        v = v - (h) * ((k * x / m) + (k * xint / m)) / 2.0;
        x = x + (h) * (v + vint) / 2.0;
        E = (m * v * v + k * x * x) / 2.0;
        t += h;
    }

 A: It works for me (using this method): 

This is the spreadsheet I used (energy column not shown):

The prediction is plain Euler: C2 = A2+B2*0.1, D2 = B2+0.1*(-A2).
The correction is trapezoidal: A3 = A2+0.1/2*(B2+D2), B3 = B2+0.1/2*(-A2-C2) 
And so on. The way I typed the formulas in, the prediction on line 2 refers to the same time moment as row 3, which perhaps isn't logical. It could be placed one row below. But this is a matter of spreadsheet style. The point is, the method works. 
The energy does grow slowly: at the end of the plot shown above, it is $12.6$ versus initial $12.5$. But this is perfectly normal for a numerical method. 
A: You can use the (velocity) Verlet (Newton-Stormer-...) method to obtain stable energy levels over long periods of time.
There are also higher-order geometric (or symplectic) integration methods, but they have sometimes baffling coefficients.
See the text "Moving Stars Around" or the paper by Hairer-etal. for details.
A: $\newcommand{\+}{^{\dagger}}
 \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
 \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
 \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
 \newcommand{\dd}{{\rm d}}
 \newcommand{\down}{\downarrow}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
 \newcommand{\fermi}{\,{\rm f}}
 \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
 \newcommand{\half}{{1 \over 2}}
 \newcommand{\ic}{{\rm i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\isdiv}{\,\left.\right\vert\,}
 \newcommand{\ket}[1]{\left\vert #1\right\rangle}
 \newcommand{\ol}[1]{\overline{#1}}
 \newcommand{\pars}[1]{\left(\, #1 \,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\pp}{{\cal P}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
 \newcommand{\sech}{\,{\rm sech}}
 \newcommand{\sgn}{\,{\rm sgn}}
 \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}
 \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
 \newcommand{\wt}[1]{\widetilde{#1}}$
You can write your equations as first order matrix equations:
\begin{align}
\vec{\xi}\pars{t} \equiv {{\rm x}\pars{t} \choose {\rm v}\pars{t}}\,,\qquad
\totald{\vec{\xi}\pars{t}}{t} = {\cal V}\vec{\xi}\pars{t}\,,\qquad
{\cal V} \equiv \pars{%
\begin{array}{cc}
0 & 1
\\[2mm]
-\,{k \over m} & 0
\end{array}}
\end{align}

$$
\xi_{n + 1} = {\cal R}\xi_{n}\,,\qquad
{\cal R} \equiv \pars{%
\begin{array}{cc}
1 & \Delta t
\\[2mm]
-\,{k \over m}\,\Delta t & 1
\end{array}}
\quad\mbox{with}\quad\left\vert%
\begin{array}{rcl}
\quad\vec{\xi}_{n} & \equiv & \vec{\xi}\pars{t}\,,\quad t \equiv n\Delta t
\\[2mm]
\quad\vec{\xi}_{n + 1} & \equiv & \vec{\xi}\pars{t + \Delta t}\,,\quad\Delta t > 0
\end{array}\right.
$$
  and $\ds{n = 0,1,2,\ldots}$. In order to guarantees stability we need
  $\ds{\Delta t < \root{m \over k}}$.

