clc; clear all; close all;
% Reactor batch donde se da  A --> B --> C

k1 = 1; k2 = 0.5; % min-1

function y = xprima(t,x,par)
  Ca = x(1); Cb = x(2); Cc = x(3); % para facilitar reconocer las variables
  k1=par(1); k2=par(2); % para facilitar reconocer los parámetros
  dCadt = -k1*Ca;
  dCbdt = k1*Ca - k2*Cb;
  dCcdt = k2*Cb;
  y = [dCadt dCbdt dCcdt];
endfunction

Cao = input('Cao = ') ;% mol/L
Cbo = 0; Cco = 0;
xo = [Cao Cbo Cco]';
t = 0:0.1:10;
par=[k1 k2];
x = lsode(@(x,t) xprima(t,x,par), xo,t);

plot(t,x)
xlabel('tiempo (min)')
ylabel('conc. (mol/L)')
legend('Ca','Cb','Cc')