function MultipleDEQexample
clear, clc, format short g, format compact
tspan = [0 3.]; % Range for the independent variable 
y0 = [1.; 0; 0]; % Initial values for the dependent variables 
disp(' Variable values at the initial point ');
disp([' t    = ' num2str(tspan(1))]);
disp('           y                  dy/dt         ');
disp([y0 ODEfun(tspan(1),y0)]);
[t,y]=ode45(@ODEfun,tspan,y0);
for i=1:size(y,2)
    disp([' Solution for dependent variable y' int2str(i)]);
    disp(['              t                  y' int2str(i)]);
    disp([t y(:,i)]);
    plot(t,y(:,i));
    title([' Plot of dependent variable y' int2str(i)]);
    xlabel(' Independent variable (t)');
    ylabel([' Dependent variable y' int2str(i)]);
    pause
end
%- - - - - - - - - - - - - - - - - - - - - -
function dYfuncvecdt = ODEfun(t,Yfuncvec); 
A = Yfuncvec(1); 
B = Yfuncvec(2); 
C = Yfuncvec(3); 
k1 = 1; 
k2 = 2; 
dAdt = -k1 * A; 
dBdt = k1 * A - k2 * B; 
dCdt = k2 * B; 
dYfuncvecdt = [dAdt; dBdt; dCdt];