multi variable function for a vector

Hi. I want to solve the set of following equations:

V=pi*r^2*L
I=K*r^2/L^2*((L^2+r^2)^0.5)-r

in which r& L are the answers. I have to solve the problem for a set of V as volume and find the graph for it.

I have written the following function which uses a gauss elimination in a newton raphson approach to find the answer. The function works properly for just one V. yet I want to change it so that I can find a matrix [x,2] for a range of v=0.001:0.001:0.01

I would appreciate any urgent help Thanks

%%N.R
function [x,ea, iter] = Hamedani_NewtRaph(func,vars,x0,es,maxit)
%function [ x, ...] = Hamedani_NewtRaph(func,var,x0,es,maxit)
%delivers the value of x, ea, iter for a system of non-linear equations based on
%Newton Raphson method while input functions are symbolic functions.
% INPUTS:
% func: symbolic vector of our functions
% vars: symbolic variables used in func
% x0= initial guess for our solution
% es= Acceptable relative error (if not set, the default is 0.0001)
% Maxit= Maximum number of iteration (if not set, the default is 50)
% OUTPUTS:
% x: the vector which represents the final length and radius
% ea= maximal approximate error in final solution
% iter: the number of iterations to reach the acceptable solution
if nargin<2, error ('at least 2 arguments must be entered'); end
if nargin<3 || isempty(es), es=0.0001; end
if nargin<4 || isempty(maxit), maxit=50; end
% primary values:
iter = 0; ea=1; x=x0; 
while (1)
  % to find the jacobian
  J = jacobian(func, vars);
  % to assign the values to the jacobian
  Jval = eval (subs(J, vars, x));
  f = eval (subs(func, vars, x));
  % for the iterations, rather than  dx=Jval\f we can call the Naive Gauss function
  dx = Hamedani_Gauss(Jval,-f)';
  x = x+dx;
  iter = iter+1;
  ea = max(abs(dx./x));
  if iter>maxit|| ea<=es, break, end
end
end