So you meant to solve dT/dt = 0.25(75-T) with initial condition of T(0) = 0? If so, then the exaction solution is T(t) = 75-75e^(-0.25t). Let me know if this is what you wanted.
function [y, t,N] = rk4(f,to,tf,ya,h)
N= ceil((tf-to)/h);
halfh = h / 2;
y(1,:) = ya;
t(1) = to;
h6 = h/6;
for i = 1 : N
t(i+1) = t(i) + h;
th2 = t(i) + halfh;
s1 = f(t(i), y(i,:));
s2 = f(th2, y(i,:) + halfh * s1);
s3 = f(th2, y(i,:) + halfh * s2);
s4 = f(t(i+1), y(i,:) + h * s3);
y(i+1,:) = y(i,:) + (s1 + s2+s2 + s3+s3 + s4) * h6;
end
%% This is your mfile. Run this alone after you have saved the above function in the same
% folder as the mfile.
yo = 0;
h = 0.1;
to = 0;
tf = 24;
Exact =@(t)(75-75*exp(-.25*t));
f = @(t,y)(0.25*(75-y));
[y, t,N] = rk4(f,to,tf,yo,h);
plot(t,Exact(t),'linewidth',2.5)
hold on
plot(t,y,'--','linewidth',2.5)
b = xlabel('t');
set(b,'Fontsize',15);
b = ylabel('y');
set(b,'Fontsize',15);
a= title('Exact Solution and Numerical solution');
set(b,'Fontsize',15);
legend('Exact Soln.','Numerical Soln.','location','best')
grid
