implementing Euler method in matlab for second order ODE

Question

I have to use the Euler method for the differential equation : $$\begin{cases} x^{\prime}=y \\ y^{\prime}=-\frac{k}{m}x-\frac{\beta}{m}x^{3} \end{cases}$$ with $k=4, \beta =-0.04 , m=1$ in matlab. We already got the code:

h=0.1; 
Tmax=50;
n=Tmax/h;
t=[0:h:Tmax];
x(1)=a;
y(1)=b;
for i=1:n
    x(i+1)=x(i)+h*f(t(i),x(i),y(i));
    y(i+1)=y(i)+h*g(t(i),x(i),y(i));
end
plot(t,x)
plot(x,y)

I replaced f and g by : y(i) and -4x(i)+0.04x(i)^3 but this did not work.

---

(error from comment) we have to choose some starting values and stepsize. When i choose $h=0.1$, $x(1)=5$, $y(1)=0$ i get the errors :

Error using plot 
    Vectors must be the same length. Error in Untitled (line 13) plot(t,x). 

Is there a reason for that? or are there any interesting of not usefull starting points, stepsizes?

---

Could someone help me?

Answer 1

The problem is that floating point arithmetic is not exact. Thus the array t=[0:h:Tmax] does not necessarily contain n+1 elements where n=Tmax/h, as you assume for the first plot. Print out the lengths of t,x and t(end) to see this directly.

Due to the floating point noise it may happen that n*h>Tmax so that the last element of t is around Tmax -h. $h=0.1$ has a binary representation that evaluates to $0.10000000000000000555$, thus $500$ times adding $h$ is computed as $50.00000000000044053649$, which is larger than $T_{max}=50$.

The safest way is to measure the length of t directly, n=length(t)-1. To ensure that Tmax is reached, the time can be constructed as t=[0:h:Tmax+0.99*h].

Or compute n=Tmax/h, apply proper rounding to get an integer, and set t=[0:n]*h.