MATLAB: 4th order runge kutta for spring mass sytem

runge kuttaspring mass systemvibrations

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)

Best Answer

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));
Related Question