4th order runge kutta for spring mass sytem

What is wrong with the code: (Also, I am a beginner, so any suggestions for where can i learn matlab for solving odes and pdes without using ode45)

function rk()
Fo=1;m=1;wn=1;w=0.4;z=0.03;
f1=@(t,y1,y2)(y2);
f2=@(t,y1,y2)(Fo/m*sin(w*t)-2*z*wn*y2-wn*y1^2);
h=0.01;
t(1)=0;y1(1)=0;y2(1)=0;
for i=1:10000
    t(i+1)=t(i)+h;
    k1y1 = h*f1(t(i),    y1(i),       y2(i));
    k1y2 = h*f2(t(i),    y1(i),       y2(i));
    k2y1 = h*f1(t(i)+h/2,y1(i)+k1y1/2,y2(i)+k1y2/2);
    k2y2 = h*f2(t(i)+h/2,y1(i)+k1y1/2,y2(i)+k1y2/2);
    k3y1 = h*f1(t(i)+h/2,y1(i)+k2y1/2,y2(i)+k2y2/2);
    k3y2 = h*f2(t(i)+h/2,y1(i)+k2y1/2,y2(i)+k2y2/2);
    k4y1 = h*f1(t(i)+h,  y1(i)+k3y1,  y2(i)+k3y2);
    k4y2 = h*f2(t(i)+h,  y1(i)+k3y1,  y2(i)+k3y2);
    y1(i+1)=y1(i) + (k1y1 + 2*k2y1 + 2*k3y1 + k4y1)/6;
    y2(i+1)=y2(i) + (k1y2 + 2*k2y2 + 2*k3y2 + k4y2)/6;
end
plot(t,y1)

[Sugar Daddy]

[darova]

everything is ok, just redure iterations steps

for i=1:1000

[madhan ravi]

I suggest you to watch Cleve Molers Videos about odes.

[Dhrumil Patadia]

thank you

[Dhrumil Patadia]

Why is the number of iteration steps causing a problem. And what if i want solution for a longer time period?

Answer 1

[Alan Stevens]

I think your definition of f2 is causingthe problem. Shouldn't it be:

f2=@(t,y1,y2)(Fo/m*sin(w*t)-2*z*wn*y2-wn*y1*abs(y1));

[Dhrumil Patadia]

i don't understand the difference between y1^2 & y1*abs(y1)

That does solve the problem though!

thanks

[Alan Stevens]

When y1 is positive there is no difference; however, when y1 is negative, y1^2 is posiive, but y1*abs(y1) is negative.

[Dhrumil Patadia]

Oh!!

its supposed to be wn^2 * y1, i squared the wrong term. sorry

thanks for your help anyways