FUNCTION FILE
function dqdt=divinefun(t,q)
dqdt(1)= (-15/500*q(1)) + (5/1000*q(2)) + 1;
dqdt(2)= (15/500*q(1)) - (18/1000)*q(2) + (3/400*q(3));
dqdt(3)= (13/1000*q(2)) - (13/400*q(3));
dqdt= dqdt'

SIMULATION FILE
commandwindow
clear
clc
close all
[t,q]=ode45('divineodefunwith3eqn',[0:40:1200],[0 0 0]);
table=table(t,q)

figure(1)
subplot(3,1,1)
plot(t,q(:,1),'green-o')
xlabel('Time(min)')
ylabel('Volume(Litre)')
legend('Tank1')
grid on
subplot(3,1,2)
plot(t,q(:,2),'blue-.*')
xlabel('Time(min)')
ylabel('Volume(Litre)')
legend('Tank2')
grid on
subplot(3,1,3)
plot(t,q(:,3),'red-.+')
xlabel('Time(min)')
ylabel('Volume(Litre)')
legend('Tank3')
grid on