%% p003 DiffFinMCK
%prog001: Solucion numerica Diferencias Finitas
%mx''+ c*x'+ k*x = f where x(0) = 0.01 and x'(0) = 0

clear all;
close all;
clc;

%% datos iniciales

m = 1;          %masa (kg)
c = 0.5;        %amortiguador viscoso (kg/s)
k = 40;        %resorte (N/m)
ti = 0;         %tiempo incial (s)
tf = 100;       %tiempo final  (s)
dt = 0.01;      %intervalo de tiempo (s)
t = ti:dt:tf;   %vector de tiempo (s)
xmin = -2.5;    %minimo
xmax = 2.5;     %maximo

x1 = 0;         %posicion incial (m)
v1 = 0;         %velocidad inicial (m/s)

%dimension tiempo
Nt = length(t);
N = Nt-1; 

%fuerza
tf1 = 10;
tf2 = 30;
Nf1 = round(tf1/dt);
Nf2 = round(tf2/dt);

f(1:Nt) = 0;            %fuerza inicializada con ceros
f(Nf1:Nf2) = 50;        %fuerza incializada en parte con 50 (N)


%% proceso principal

%inicializacion de la solucion
xn = zeros(1,Nt);   

%posicion inicial t = 0
xn(1) = x1;    

%fuerza inicial t = 0
f1 = f(1);    

%aceleracion inicial t = 0
a1 = ( f1 - c*v1 - k*x1 )/m;

%posicion inicial t = dt
xn(2) = x1 + v1*dt + 1/2*a1*dt^2;

%coeficientes 
cof1 = m/dt/dt + c/2/dt;
cof2 = 2*m/dt/dt - k; 
cof3 = m/dt/dt - c/2/dt;

%loop para los demas puntos
for j = 2:N 

    xn(j+1) = ( f(j) + cof2*xn(j) - cof3*xn(j-1) )/cof1;
        
end;

%% resultados
      
figure(1)
plot(t,xn);
title('solucion de diferencias finitas');
xlabel('tiempo t (s)');
ylabel('posicion x(t) (m)');
grid on;
box on;