The first order ODE is defined below

% Example of a first order system
% P. Venkataraman
% dht/dt + h/(AR) = Qin(t)/A; h(0) = 40
%
% A = 78.5, R = 100;
% (a) Qin(t) = 0; (natural response)
% (b) Qin(t) = 0.3(step response)
%-----------------------------------
% Natural response and forced response
%--------------------
clear      % clears all variables
clc      % clears window and starts at top
A = 78.5; R = 100;  Qina = 0.3;

% dsolve is used to solve differential equations
% by default independent variable is t ot time
h = dsolve('Dh = -h/(A*R)','h(0)=40');  % natural solution
ha = dsolve('Dh = -h/(A*R)+ Qina/A','h(0)=40');  % forced solution (step)

hn = subs(h);  % substitutes available parameter values in solution
fprintf('\nNatural Motion -Case (a) ')
% natural motion solution printed iin Command Window
pretty(hn)
fprintf('\n\nForced Motion - Step - Case (b)')
haa = subs(ha);
% step response printed to Command window
pretty(haa)

format compact
set(gcf,'Menubar','none','Name','First Order ODE - Symbolic Solution', ...
    'NumberTitle','off','Position',[10,350,800,350], ...
    'Color',[0.5 0.7 0.7]);
tr = 0:100:5000;
for i = 1:length(tr)
    t = tr(i);  % by default t is the independent variable
    han(i) = subs(h);
    hap(i) = subs(ha);
end
%-----------------------------------------
hp(1) = subplot('Position',[0.1 0.15 0.35 0.75]);
plot(tr,han,'r')
title('\bf\itNatural Motion - No Input','Color','k', ...
    'VerticalAlignment','bottom')
ylabel('h(m)', ...
    'FontWeight','b','Color','b', ...
    'VerticalAlignment','bottom');

xlabel('time(s)', ...
    'FontWeight','b','Color','b', ...
    'VerticalAlignment','top');  % boldface blue
legend('Natural Motion');
grid

%-----------------------------------------
hp(2) = subplot('Position',[0.55 0.15 0.35 0.75]);
plot(tr,hap,'k--')
title('\bf\itForced Motion - Constant Input','Color','k', ...
    'VerticalAlignment','bottom')
ylabel('h(m)', ...
    'FontWeight','b','Color','b', ...
    'VerticalAlignment','bottom');

xlabel('time(s)', ...
    'FontWeight','b','Color','b', ...
    'VerticalAlignment','top');  % boldface blue
legend('Forced Response');
grid
txtstr = strcat('Inflow Q_{in}(t)[m^3/s] =  ',num2str(Qina));
text(tr(5),min(hap)+1,txtstr);