finding terminal velocity using newton raphson

We are trying to find the terminal velocity of a sphere falling in water

g= 9.8; % in m/s^2
rau_p=  7850; %particle density in Kg/m^3
rau_f= 1000; % fluid density in Kg/m^3
vis_f= .001; % fluid viscoisty in Pa.s
C_D=zeros(1); %initial guess
plotdata_x=[]; 
plotdata_y=[];
for Dp=0.0001:.00001:0.1 %in m
    Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f)); %terminal velocity in m/s
    Re=rau_f*Vt*Dp/vis_f;
    if Re<=.1 
        C_D=24/Re; 
        Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
    elseif Re>0.1 && Re<=1000
        C_D= (24/Re)*(1+.14*(Re^.7));
        Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
    elseif Re>1000 && Re<=350000
        C_D=0.445;
        Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
    elseif Re>350000 && Re<=1000000
        C_D=interp1([350000 400000 500000 700000 1000000],[.396 .0891 .0799 .0945 .11],Re);
        Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
    else %Re>1000000
        C_D= .19-((8*10^4)/Re);
        Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
       
    end
     plotdata_x(end+1)=Dp;
     plotdata_y(end+1)=Vt;
       
end
plot(plotdata_x,plotdata_y) 
xlabel('D(m)')
    ylabel('Vt(m/s)')

This method is working fine however, we are also required to find the terminal velocity using Newton-Raphson method and this is what we have done so far but it is not working

This is the newton-raphson code

    function Xs = Newton(f,Xest,Err,imax)
%Newton finds the root of f=0 near the point X0 using Newton's method

%Input variables:
% f     : Name of a user-defined function that calculates function for a given x
% Xest  : Initial estimate of the solution
% Err   : Relative approximate error
% imax  : Maximum number of iterations

%Output variables: 
% Xs    : Solution
syms x;
fprintf('The derivative function is:')
g = diff(f) %The Derivative of the Function

disp('Iter    Xi       Relative approximate error');
for i = 1:imax
    Xi = Xest - vpa(subs(f,x,Xest))/vpa(subs(g,x,Xest));
    fprintf('%2i \t %f \t %f \n', i, Xi, abs((Xi - Xest)/Xest));
    if abs((Xi - Xest)/Xest) < Err
        Xs = Xi;
        break
    end
    Xest = Xi;
end
if i == imax
    fprintf('Solution was not obtained in %i iterations.\n',imax)
    Xs=('No answer');
end

and this is what we tried

clc
clear
Dw=1000;     %Density of water in kg/m^3
Dp=7850;    %Density of particle(iron) in kg/m^3
U=0.001;     %Viscosity in kg/m.sec or pa.sec
g=9.807;    %Accelaration due to gravity in m/sec^2
vt=1;        %Initial guess in m/s
dp=input('Enter dp Value: ');
Rep=(dp*vt*Dw)/U;%Reynolds number of particle
for i=1:1000
    if Rep<0.1
        Rep=(Dw*vt*dp)/U;
        cd=24/Rep;
        vtd=6713/360000;        %Derivative of the Vt in laminar region
        V=((4*(Dp-Dw)*g*Dw)/(3*cd*Dw))^0.5;
        vt=vt-(V/vtd);          %Newton Rapshon Method
    elseif (0.1<Rep)&&(Rep<=800)
        Rep=(Dw*vt*dp)/U;
        cd=(24*(1+0.15*Rep^0.687))/Rep;
        vtd=6713/360000;        %Derivative of the Vt in laminar region
        V=((4*(Dp-Dw)*g*Dw)/(3*cd*Dw))^0.5;
        vt=vt-(V/vtd);          %Newton Rapshon Method
    elseif (1000<Rep)&&(Rep<=350000)
        Rep=(Dw*vt*dp)/U;
        cd=0.445;
        V=((4*(Dp-Dw)*g*Dw)/(3*cd*Dw))^0.5;
        vt=V;
    else
        Rep=(Dw*vt*dp)/U;
        cd=0.19-((8*10^4)/Rep);
        V=((4*(Dp-Dw)*g*Dw)/(3*cd*Dw))^0.5;
        vtd=-161112/(5*vt^2*(2400/vt - 57/100)^2);
        vt=vt-(V/vtd);
    end
end

[Dhananjay Kumar]

Can you please specify what you mean by 'not working' ? Also please write your purpose in terms of mathematical equations and operations so that it's easy to understand for anyone.