syms T(t) C(t)
con1=T(0)==16 ;
con2=C(0)==1 ;
ode1=diff(C) == -exp(-10./(T+273)).* C;
ode2=diff(T) == 1000*exp(-10./(T+273)).* C-10*(T-20);
vars = [T(t) ; C(t)]
V = odeToVectorField([ode1,ode2])
M = matlabFunction(V,'vars', {'t','Y'})
interval = [0 5];
y0 = [16 1];
ySol = ode45(M,interval,y0);
tValues = linspace(interval(1),interval(2),1000);
yValues = deval(ySol,tValues,1);
plot(tValues,yValues,'-o')
figure
yValues = deval(ySol,tValues,2);
plot(tValues,yValues,'-or')