# implementing Euler method in matlab for second order ODE

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?

• What exactly did not work? I'd suspect that the array lengths are not exactly equal and that's why the first plot command fails. Change to t=[0:h:Tmax]; n = len(t)-1; to get equal lenghts. You might also investigate the plotting functions to make sure that the phase plot is in a different canvas/figure than the component plot. That the multiplication should be written as -4*x(i)+0.04*x(i)^3 should not have to be mentioned. – LutzL Nov 11 '18 at 9:56
• h=0.1; Tmax=50; n=Tmax/h; t=[0:h:Tmax]; n = len(t)-1; x(1)=1; y(1)=2; for i=1:n x(i+1)=x(i)+h*y(i); y(i+1)=y(i)+h*(-4*x(i)+0.04*x(i)^3); end plot(t,x) plot(x,y) this is what i have now. I get the error : "Undefined function 'len' for input arguments of type 'double'. Error in Untitled (line 7) n = len(t)-1;" – B.D Nov 11 '18 at 10:01
• Then use length to get the length of the array. And you only need to compute the final value of n, the first one is not used. Please add the original error (or its description if no error message) to the question text. As a coding problem, the better place would be the stackoverflow.com forum. – LutzL Nov 11 '18 at 10:06
• 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 usfull starting points,stepsizes? Thanks for helping – B.D Nov 11 '18 at 10:26

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.