You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: flex_analytical.m
    path of file: ./repos/gmtmex/test/potential
    
    the code of the file until where you have to start completion: function flex_analytical
% Script to calculate the Matlab exact solutions to disc flexure so that we
% can compare it to gravfft 3-D crossections.  The disc dimensions are
h = 4000; r = 15; rhol = 2800; rhoi = 2400; rhom = 3300; rhow = 1000;
Te = 10;
x = (0:512)';
% a) Regular single density case
wa = flexdisc (x, h, r, rhol, rhom, rhol, rhow, Te);

    
function flex_analytical
% Script to calculate the Matlab exact solutions to disc flexure so that we
% can compare it to gravfft 3-D crossections.  The disc dimensions are
h = 4000; r = 15; rhol = 2800; rhoi = 2400; rhom = 3300; rhow = 1000;
Te = 10;
x = (0:512)';
% a) Regular single density case
wa = flexdisc (x, h, r, rhol, rhom, rhol, rhow, Te);
% b) Exact dual density case
wb = flexdisc (x, h, r, rhol, rhom, rhoi, rhow, Te);
% c) Approximate single density case
A = sqrt ((rhom - rhoi)/(rhom-rhol));
wc = A * flexdisc (x, h, r, rhoi, rhom, rhoi, rhow, Te);
figure(1);  clf;
plot (x, wb, 'g.', 'MarkerSize', 10); hold on
plot (x, wa, 'b-');
plot (x, wc, 'r-');
tmp = [x wa wb wc];
save flex_analytical.txt tmp -ascii -tabs

function w = flexdisc (x, h, r, rhol, rhom, rhoi, rhow, te)
% FLEXDISC calculates the flexural deformation caused by a
% disc of given dimensions. Infill density may differ from
% load density
%
% w = flexdisc (x, h, r, rhol, rhom, rhoi, rhow, te)
%
%     x		- where to evaluate flexure (all x >= 0) (km)
%     h		- height of disc (m)
%     r		- radius of disc (km)
%     rhol	- density of load
%     rhom	- density of mantle
%     rhoi	- density of infill
%     rhow	- density of water
%     te	- elastic plate thickness (km)
%
% On output, w is the deflections.

te = 1000.0 * te;
yi = flex3dk (te, rhom - rhol);		% Get inside flex wavenumber
yo = flex3dk (te, rhom - rhoi);		% Get outside flex wavenumber
airy = (rhol - rhow) / (rhom - rhol);	% Get density balance
in = find (x <= r);			% Get indices to points inside/outside
out = find (x > r);
xm = 1000.0 * x;			% Convert to meters
x(in) = xm(in) * yi;			% Scale to flexural wavelengths here instead of later
x(out) = xm(out) * yo;
ri = 1000.0 * r * yi;
if max(x(:)) > 8.0
	disp(['Warning: Solution inaccurate for distances beyond r = ' num2str(8/(1000.0*yo)) ' km'])
end