Matlab

a=0;
t=linspace(a,b,N);
dt=(b-a)/N;
w(1)=0;
teta(1)=1;
X(1)=10;
vX(1)=0;
for i=1:N-1alfa(i)=(-(m+M)*g*sin(teta(i))+k*X(i)*cos(teta(i))-m*l*cos(teta(i))*sin(teta(i))*w(i)^2)/(l*(m+M-m*(cos(teta(i)))^2));w(i+1)=w(i)+dt*(-(m+M)*g*sin(teta(i))+k*X(i)*cos(teta(i))-m*l*cos(teta(i))*sin(teta(i))*w(i)^2)/(l*(m+M-m*(cos(teta(i)))^2));%w=velocidad en tetavX(i+1)=vX(i)+dt*((-k*X(i)+m*(l*sin(teta(i))*w(i)^2-l*cos(teta(i))*alfa(i))/(M+m)));
teta(i+1)=teta(i)+dt*w(i);
X(i+1)=X(i)+dt*vX(i);
endx=X+l.*sin(teta);
y=-l.*cos(teta);
figure(1)
subplot(2,1,1)
plot(t,X)
xlabel ('t(s)')
ylabel('X(m)')
subplot(2,1,2)
plot(t,teta)xlabel('t(s)')
ylabel('teta(rad)')
figure(2)
plot(x,y)
xlabel ('x(m)')
ylabel ('y(m)')
figure (3)
subplot(2,1,1)
plot(t,vX)
xlabel('t(s)')
ylabel ('vX (m/s)')
subplot(2,1,2)
plot(t,w)
xlabel ('t(s)')
ylabel ('w (rad/s)')

%Ahora hallo la energía para ver si seconserva:
figure(6)
E=+(1/2).*k.*X.^2+m.*g.*y+(1/2).*M.*(vX.^2)+(1/2).*m.*(vX.^2)+m.*vX.*l.*cos(teta).*w+(1/2).*m.*(l.^2).*(w.^2);plot(t,E./E(1)-1)
xlabel ('t(s)')
ylabel ('E/E(1)-1')

figure(7)plot(t,(1/2).*m.*(l.^2).*(w.^2),t,(1/2).*k.*X.^2,t,+m.*g.*y,t,(1/2).*M.*(vX.^2),'k',t,(1/2).*m.*(vX.^2),t,m.*vX.*l.*cos(teta).*w)
legend('1', 'EP M', 'EP Y', 'EC M', '5', '6')
%%
%Runge-Kutta