EDITED
If you want to pass a new parameter to a function parametrize your function
% Set extra parameters
tf = 72; % h
h0 = 2; % m
g = 9.81;
a=input('leakage value : ');
L=a*g*h;
%prompt = 'What is the leakage value? ';
%a = input(prompt);
%L = a*g*h;
% Solve the DAE
opts = odeset('Reltol',1e-5);
[t,h] = ode45(@(t,h)tankf1(t,h,L),[0 tf],h0,opts);
% Plot the graph
plot(t,h)
title('Tank model solved by Method 1')
xlabel('Time (h)')
ylabel('Level (m)')
grid on
function dhdt = tankf1(t,h,L)
% TANKF1: Tank level model - Derivative function for solution method 1
% (substitution of algebraic variables)
% Tank model parameters and constants
A = 12.5; % m^2
CV1 = 3.41; % (m^3/h)/kPa^0.5
CV2 = 3.41; % (m^3/h)/kPa^0.5
P0 = 100; % kPa
P3 = 100; % kPa
rho = 1000; % kg/m^3
g = 9.81; % m/s^2
P1 = P1function(t); % kPa
% Algebraic equations for F1, F2 and P2 have been substituted into the
% differential equation for dh/dt to create a pure ODE system:
dhdt = 1/A*(CV1*sqrt(P1-P0-rho*g*h/1000) - CV2*sqrt(P0+rho*g*h/1000-P3))-L;
end
function P1 = P1function(t) % should be defined at last
% P1FUNCTION: Upstream pressure (kPa) as a function of time (h)
if t<=10
P1 = 139.5;
elseif t<=11
P1 = 139.5 + 60.5*(t-10);
elseif t<=40
P1 = 200;
elseif t<=55
P1 = 200 - 30*(t-40)/15;
else % t>55
P1 = 170;
end
end
example of the graph produced:
