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);
    z(i+1)=z(i)+h*((-1.24*10^(-4))*z(i));
    t(i+1)=t(i)+h;
    i=i+1;
end
[t',M',z']
hold on
plot(t,M)
plot(t,z,'--r')
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)
%Método de Euler
h=0.01;
M0=1;
N0=100;
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
M(1)=M0;
N(1)=N0;
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
hold on
plot(t2,N,'--g')
plot(t1,M)
hold off
disp(t1(end))
disp(t2(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).

%Trapecio
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;
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)

%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