Terceira

clear;
constante = 50;
m = 4.5;
k = 1400;
f1 = @(t, y1, y2) (y2);
f2 = @(t, y1, y2) -( (constante/m) * y2 + (k/m) * y1);
t0 = 0;
tN = 5;
y1_0 = 0.02;
y2_0 = 0;
max_it = 200;
h = (tN - t0)/max_it;
t = linspace(t0, max_it,max_it+1);
w1 = linspace(0, max_it ,max_it+1);
w2 = linspace(0, max_it ,max_it+1);
y1 = linspace(0, max_it ,max_it+1);
y2 = linspace(0, max_it, max_it+1);
t(1) = t0;
y1(1) = y1_0;
y2(1) = y2_0;
w1(1) = y1_0;
w2(1) = y2_0;
printf("\t\t\t\t w1\t\t\t w2\t\t\t \n\n");
for i = 1:max_it
k1 = h * f1( t(i), w1(i), w2(i) );
l1 = h * f2( t(i), w1(i), w2(i) );
k2 = h * f1( t(i) + h/2 , w1(i) + k1/2 , w2(i) + l1/2);
l2 = h * f2( t(i) + h/2 , w1(i) + k1/2 , w2(i) + k1/2);
k3 = h *f1 ( t(i) + h/2, w1(i) + k2/2, w2(i) +l2/2 );
l3 = h *f2 ( t(i) + h/2, w1(i) + k2/2, w2(i) +l2/2 );
k4 = h * f1( t(i) + h, w1(i) + k3, w2(i) + l3/2 );
l4 = h * f2( t(i) + h, w1(i) + k3, w2(i) + l3/2 );
w1(i+1) = w1(i) + (k1 + 2 *(k2 + k3) + k4)/6;
w2(i+1) = w2(i) + (l1 + 2 *(l2 + l3) + l4)/6;
printf("%.1f\t\t %.7f\t\t %.7f\t\t\n ",t(i+1),w1(i+1), w2(i+1));
end
plot(t, w1, 'b', t , w2, 'r');
legend('Posicao','Velocidade');