clear all, clf
h=input(' Introduce el valor del paso: ')
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t(1)=t0;
M(1)=M0;
z(1)=M0;
i=1;
j=1;
k=1.24*10^-4;
while M(i)>0.08*M0;
    M(i+1)=((1-k*h/2)*M(i))/(1+k*h/2);%Trapecio
    z(i+1)=z(i)+h*((-1.24*10^(-4))*z(i));%Euler
    t(i+1)=t(i)+h;
    i=i+1;
end
[t',M',z']
hold on
plot(t,M,'--r')
plot(t,z)
hold off
disp(t(end))%Muestra en pantalla el valor del vector "t" de tiempo

clear all, clf
%M'(t)=-1,24*10^-4*M(t)
% A continuación exponemos una comparación de las dos gráficas que resultan
% de variar el paso de 0.1 a 0.01
h1=0.1;
h2=0.01;
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t1(1)=t0;
t2(1)=t0
M1(1)=M0;
M2(1)=M0;
i=1;
n=1;
%El primer Bucle utiliza un paso de 0.1
while M1(i)>(0.08*M0);
    M1(i+1)=M1(i)+h1*((-1.24*10^(-4))*M1(i));
    t1(i+1)=t1(i)+h1;
    i=i+1;
end
%El segundo bucle utiliza un paso de 0.01
while M2(n)>(0.08*M0);
    M2(n+1)=M2(n)+h2*((-1.24*10^(-4))*M2(n));
    t2(n+1)=t2(n)+h2;
    n=n+1;
end
%expone en pantalla la antigüedad de los restos arqueológicos
hold on
plot(t1,M1,'r')
plot(t2,M2,'--b')
hold off

clear all, clf
%M'(t)=-1,24*10^-4*M(t)
%Método de Euler
h=0.01;
M0=1;
N0=100;
Z0=50;
t0=0;
t1(1)=t0;%Utilizamos el vector t1 para calcular la edad de un resto arqueológico con una supuesta masa de 1
t2(1)=t0;%Utilizamos el vector t1 para calcular la edad de un resto arqueológico con una supuesta masa de 100
t3(1)=t0;
M(1)=M0;
N(1)=N0;
Z(1)=Z0;
i=1;
%utilizamos un bucle "while" para evitar que la función siga adoptando
%valores una vez superado llegue a ese 8% que se encuentra en el hueso
while M(i)>(0.08*M0);
    M(i+1)=M(i)+h*((-1.24*10^(-4))*M(i));
    t1(i+1)=t1(i)+h;
    i=i+1;
end
j=1;
while N(j)>(0.08*N0);
    N(j+1)=N(j)+h*((-1.24*10^(-4))*N(j));
    t2(j+1)=t2(j)+h;
    j=j+1;
end
n=1
while Z(n)>(0.08*Z0);
    Z(n+1)=Z(n)+h*((-1.24*10^(-4))*Z(n));
    t3(n+1)=t3(n)+h;
    n=n+1;
end
hold on
plot(t2,N,'g')
plot(t1,M)
plot(t3,Z,'r')
hold off
disp(t1(end))
disp(t2(end))
disp(t3(end))
%Sacamos en pantalla el valor de t1 y t2 para comparar la edad de ambos.
Error=t1-t2;
disp(Error(end))%Expone en pantalla la diferencia entre t1 y t2 y queda 
%demostrado que no depende de la cantidad inicial de M (materia).

%El método del trapecio es 'indirecto' lo que quiere decir que previamente
%debemos despejar analíticamente el valor de f(t(n+1),y(n+1)
clear all, clf
h=input(' Introduce el valor del paso: ')
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t(1)=t0;
M(1)=M0;
i=1;
k=1.24*10^-4;
%Igual que en el caso anterior, Utilizaremos un bucle "while" para que el
%programa se detenga automáticamente una vez reducida la cantidad de C14
while M(i)>0.08*M0;
    M(i+1)=((1-k*h/2)*M(i))/(1+k*h/2);
    t(i+1)=t(i)+h;
    i=i+1;
end
disp(t(end))
plot(t,M)

clear all, clf
h=input(' Introduce el valor del paso: ')
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t(1)=t0;
M(1)=M0;
z(1)=M0;
i=1;
j=1;
k=1.24*10^-4;
%La "M" representa los valores obtenidos por el método de Euler, y la "z" el
%vector de volores obtenidos mediante Trapecio
while M(i)>0.08*M0;
    %Trapecio
    M(i+1)=((1-k*h/2)*M(i))/(1+k*h/2);
    %Euler
    z(i+1)=z(i)+h*((-1.24*10^(-4))*z(i));
    t(i+1)=t(i)+h;
    i=i+1;
end
%A continuación se muestra una tabla comparativa para los distintos valores
%de "t" comparando los obtenidos mediante "Euler" y mediante "Trapecio"
[t',M',z'];
hold on
plot(t,M,'--r','linewidth',1.5)
plot(t,z,'linewidth',1.5)
hold off

%Rounge Kutta
clear all, clf
%M'(t)=-1,24*10^-4*M(t)
h=input(' Introduce el valor del paso: ')
M0=input(' Introduce la cantidad inicial de Carbono 1: ')
t0=0;
t(1)=t0;
M(1)=M0;
i=1;
k=1.24*10^-4;
%Los métodos de Runge-Kutta (RK) son un conjuntos de métodos iterativos
%para la aproximación de ecuaciones diferenciales ordinarias.
while M(i)>0.5*M0;
    %K1=f(t(n),M(n))
    K1=-k*M(i);
    %K2=f(t(n)+1/2*h, M(n)+1/2K1*h)
    K2=-k*M(i)-k/2*K1*h;
    %K3=f(t(n)+1/2*h,M(n)+1/2*K2*h)
    K3=-k*M(i)-k/2*K2*h;
    %K4=f(t(n)+h,M(n)+K3*h)
    K4=-k*(M(i)+K3*h);
    %M(n+1)=M(n)+h/6*(K1+2K2+2K3+K4)
    M(i+1)=M(i)+h/6*(K1+2*K2+2*K3+K4);
    t(i+1)=t(i)+h;
    i=i+1;
end
plot(t,M,'g')
disp(t(end))

%Condiciones iniciales.
t0=0;
tf=10;
k1=5;
k2=1;
%Tamaño de paso
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%Matriz para la resolución
y=zeros(3,N+1);
%Concentraciones iniciales
A=1;B=0;C=0;
y0=[A B C]';
y(:,1)=y0;
%Resolvemos por euler partiendo de la matriz de coeficientes sabiendo que
%y'=M*y+t;
M=[-k1 0 0;k1 -k2 0;0 k2 0];
for i=1:N
    y(:,i+1)=y(:,i)+h*(M*y(:,i));
end
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
plot(t,y(3,:))
legend('Compuesto A','Compuesto B','Compuesto C','Location','best')
hold off

%Condiciones iniciales.
t0=0;
tf=10;
k1=5;
k2=1;
%Tamaño de paso
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%Matriz para la resolución
y=zeros(3,N+1);
I=eye(3);
%Concentraciones iniciales
A=1;B=0;C=0;
y0=[A B C]';
y(:,1)=y0;
%Resolvemos por euler partiendo de la matriz de coeficientes sabiendo que
%y'=M*y+t;
M=[-k1 0 0;k1 -k2 0;0 k2 0];
for i=1:N
    y(:,i+1)=(I-(h/2)*M)\(y(:,i)+(h/2)*(M*y(:,i)));%Método del trapecio
end
hold on
plot(t,y(1,:),'r')
plot(t,y(2,:),'k')
plot(t,y(3,:))
legend('Compuesto A','Compuesto B','Compuesto C','Location','best')
hold off