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.
Stability of explicit Runge-Kutta methods around a fixed point or a slow moving equilibrium is obtained when $Lh\sim 1$ with $L$ the $y$-Lipschitz constant and $h$ the step size. The exact constant on the right side depends on the order of the method, usually it is between 1 and 5. For first and second order methods that is only a guideline as the stability region does not contain a semi-circle around the origin in the negative complex half-plane.
To get a practical feeling for this, just solve your problem for various step sizes over the integration interval from $0$ to $1/3$, as $y(1/3)=0.1e^{-50}$ is a really small number, relative to the initial value and floating point precision.
def Heun(f,t0,y0,tf,h):
def HeunStep(t,y,h):
k1 = h*f(t,y);
k2 = h*f(t+h,y+k1);
return y+0.5*(k1+k2);
t, y = t0, y0
while t+1.1*h<tf: t, y = t+h, HeunStep(t,y,h)
return HeunStep(t,y,tf-t)
def f(t,y): return 49-150*(t+y)
t0,y0 = 0.0,1./3+0.1
for h in np.linspace(3.5,0.5,16)/150:
print "h=%.4e 150*h=%.4e, y(1/3)=%.12g"%(h,150*h, Heun(f,t0,y0,1./3,h))
with the results
h=2.3333e-02 150*h=3.5000e+00, y(1/3)=3382845.45109
h=2.2000e-02 150*h=3.3000e+00, y(1/3)=1820496.14682
h=2.0667e-02 150*h=3.1000e+00, y(1/3)=558708.429462
h=1.9333e-02 150*h=2.9000e+00, y(1/3)=79763.1519923
h=1.8000e-02 150*h=2.7000e+00, y(1/3)=9204.09082254
h=1.6667e-02 150*h=2.5000e+00, y(1/3)=1648.41784102
h=1.5333e-02 150*h=2.3000e+00, y(1/3)=37.6118108918
h=1.4000e-02 150*h=2.1000e+00, y(1/3)=0.740437945724
h=1.2667e-02 150*h=1.9000e+00, y(1/3)=0.00432803787737
h=1.1333e-02 150*h=1.7000e+00, y(1/3)=1.0688786734e-05
h=1.0000e-02 150*h=1.5000e+00, y(1/3)=1.14794370277e-08
h=8.6667e-03 150*h=1.3000e+00, y(1/3)=5.57823318786e-12
h=7.3333e-03 150*h=1.1000e+00, y(1/3)=2.8015784137e-15
h=6.0000e-03 150*h=9.0000e-01, y(1/3)=2.34187669257e-17
h=4.6667e-03 150*h=7.0000e-01, y(1/3)=-1.38777878078e-17
h=3.3333e-03 150*h=5.0000e-01, y(1/3)=2.60208521397e-17
As the stability region suggests, $0<h<\frac2{150}$ are valid choices in the table above, giving results close to zero. Sensible results that are practially zero as floating point numbers (and relative to the initial value) are only obtained for about $h<\frac{1.5}{150}$.
Best Answer
The general solution of this ODE is $e^{-x}+C\,e^x$. With $y(0)=1$, we conclude $C=0$. But every rounding error in any numerical method will introduce a spurious $C\neq0$, and given the notorious behavior of exponential functions, this will overwhelm your numerical approximation.