Loss of stability (unphysical energy gain) for simple pendulum equation? I am simulating a pendulum using MATLAB and noted a curious issue.
When I use zero velocity and (pi - 0.1) angular position as starting conditions for my second order ODE, the solution deviates from what physically and analytically expected after some short time (the pendulum gains energy, after Swinging back and forth it starts revolving).
The trivial code I used is as follows 
f1 = @(t,Y) [Y(2); -sin(Y(1))];
[z1,z2] = ode45(f1,[0:0.001:80],[ pi - 0.1,0])
plot(z1,z2(:,1))
Using ODE15s actually makes the phenomenon to occur earlier.
Any hint on what I am doing wrong? What is causing all this?
Thanks a lot
 A: This is not too surprising.  Energy is conserved in the true solution to the system, but the numerical approximations won't conserve it exactly.  If on the average the system gains a little bit of energy at each step, after enough steps the gain will become noticeable.
EDIT:
In this system, the energy is $y_2^2/2 - \cos(y_1)$.  The initial energy is
$-\cos(\pi - .1) \approx .995004$, very close to the energy ($1$) required to reach the vertical position, so even quite a slow accumulation of energy is enough to put you over the top.
One thing you can try, if you know a conserved quantity (energy in this case), is to use a change of variables in which that quantity is one of the new variables.  In this case with $E = y_2^2/2 - \cos(y_1)$, you know that the
solution curves should be closed curves around the origin (at which $E=-1$), 
so you can write $E + 1 = y_2^2/2 + 2\sin^2(y_1/2)$.  Then with
$\theta$ defined by $ \cos(\theta) = \sqrt{2/(E+1)} \sin(y_1/2)$ and $\sin(\theta) = y_2/\sqrt{2(E+1)}$ I get (if I haven't made any mistakes)
$$ \dot{\theta} = -\sqrt{1 - (E+1)\cos^2(\theta)/2}$$
So you solve that differential equation and take $y_1 = 2 \arcsin(\sqrt{(E+1)/2} \cos(\theta))$ and $y_2 = \sqrt{2(E+1)} \sin(\theta)$.
