function [ v ] = plot_bfun( ind )

global model = [ 	19	2.74	6.14	3.55	1000	450
			38	3	6.58	3.8	1000	450
			50	3.32	8.2	4.65	1000	60
			60	3.34	8.17	4.62	1000	60
			70	3.35	8.14	4.57	1000	80
			80	3.36	8.1	4.51	1000	100
			90	3.37	8.07	4.46	1000	100
			100	3.38	8.02	4.41	1000	100
			125	3.39	7.93	4.37	1000	150
			150	3.41	7.85	4.35	1000	150
			175	3.43	7.89	4.36	1000	150
			200	3.46	7.98	4.38	1000	150
			225	3.48	8.1	4.42	1000	150
			250	3.5	8.21	4.46	1000	150
			300	3.53	8.38	4.54	1000	150
			350	3.58	8.62	4.68	1000	150
			400	3.62	8.87	4.85	1000	180
			450	3.69	9.15	5.04	1000	180
			500	3.82	9.45	5.21	1000	250
			600	4.01	9.88	5.45	1000	450
			700	4.21	10.3	5.76	1000	500
			800	4.4	10.71	6.03	1000	600
			900	4.56	11.10	6.23	1000	1000
			1000	4.63	11.35	6.32	1000	1000
		 ];

nfreq = 30;

global Rearth = 6380.0;

global nlayers = length(model);
global nbfun = 1+(nlayers-2)*2;   # one basis function on top plus two for each level except bottom
global n2 = 2*nbfun;		# number of parameters in r1 and r2 combined

	# earth-flattening approximation

global z = -Rearth*log(ones(nlayers,1) - model(:,1)/Rearth);
global rho = model(:,2).*(ones(nlayers,1) - model(:,1)/Rearth).^5;
global vp = model(:,3).*(ones(nlayers,1) - model(:,1)/Rearth);
global vs = model(:,4).*(ones(nlayers,1) - model(:,1)/Rearth);
global Qp = model(:,5);
global Qs = model(:,6);

global myu = vs.^2 .* rho;
global lambda = vp.^2 .* rho - 2.0 * myu;
global bulk = lambda - 2.0/3.0 * myu;

global l1 = zeros(n2,1);		# current solution
global zmid = zeros(nlayers,1);
global dz = zeros(nlayers,1);

zmid(1) = 0.5*z(1);
dz(1) = z(1);

for i=2:nlayers
	zmid(i) = 0.5*(z(i)+z(i-1));
	dz(i) = z(i)-z(i-1);
endfor

z1 = [ 0:0.1:100 ];
f = z1;
g = z1;

for i=1:length(z1)
	f(i) = bfunval( ind, z1(i) );
	g(i) = bfungrad( ind, z1(i) );
endfor

figure(1);
plot(z1,f,';val;',z1,g,';grad;');

endfunction