function P2_16A
clear, clc, format short g, format compact
tspan = [0 200.]; % Range for the independent variable 
y0 = [20.; 20.; 20.]; % 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); 
T1 = Yfuncvec(1); 
T2 = Yfuncvec(2); 
T3 = Yfuncvec(3); 
W = 100; 
Cp = 2; 
T0 = 20; 
UA = 10; 
Tsteam = 250; 
M = 1000; 
dT1dt = (W * Cp * (T0 - T1) + UA * (Tsteam - T1)) / (M * Cp); 
dT2dt = (W * Cp * (T1 - T2) + UA * (Tsteam - T2)) / (M * Cp); 
dT3dt = (W * Cp * (T2 - T3) + UA * (Tsteam - T3)) / (M * Cp); 
dYfuncvecdt = [dT1dt; dT2dt; dT3dt];