% Euler Modificado
% Sistema de 3 masas y 4 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
k4=3;
m1=2;
m2=1;
m3=3;
%Condiciones iniciales
x0=-0.5;
y0=1;
z0=-0.5;
u0=1;
v0=1;
w0=1;
 
N=50; % Pasos
h=(tN-t0)/N; % Intervalo
 
 
% Primeros valores
xx(1)=x0;
yy(1)=y0;
zz(1)=z0;
uu(1)=u0;
vv(1)=v0;
ww(1)=w0;
 
% Metodo de Runge-Kutta
 
 
 
for n=1:N
    % Primera K
    k1x=uu(n);
    k1y=vv(n);
    k1z=ww(n);
    k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
    k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n)+(k3/m2)*zz(n);
    k1w=k3/m3*yy(n)-((k4+k3)/m3)*zz(n);
 
    xa=xx(n)+h/2*k1x;
    ya=yy(n)+h/2*k1y;
    za=zz(n)+h/2*k1z;
    ua=uu(n)+h/2*k1u;
    va=vv(n)+h/2*k1v;
    wa=ww(n)+h/2*k1w;
    
    % Segunda K
    k2x=ua;
    k2y=va;
    k2z=wa;
    k2u=-(k2+k1)/m1*xa+k2/m1*ya;
    k2v=k2/m2*xa-((k2+k3)/m2)*ya+(k3/m2)*za;
    k2w=k3/m3*ya-((k4+k3)/m3)*za;
 
    
    
    xx(n+1)=xx(n)+h*k2x;
    yy(n+1)=yy(n)+h*k2y;
    zz(n+1)=zz(n)+h*k2z;
    uu(n+1)=uu(n)+h*k2u;
    vv(n+1)=vv(n)+h*k2v;
    ww(n+1)=ww(n)+h*k2w;
end
hold on
t=t0:h:tN;
plot(zz+8-(xx+2),t,'-k')
hold off
axis([0,12,0,10])
hold on
title('Euler Modificado')

legend('Masa m1','Masa m2','Masa m3','Posicion de equilibrio')
hold off