Error using double (subs())

Hi,

I want to run newton's method and I am not able to convert the syms variable into double even though I am using subs.

Here is my code.

syms initial_R
u = -100;
v = @(initial_R) initial_R.*exp ((2*(1-alpha))/(2+(1-alpha)*u.*initial_R)./(2-alpha*u.*initial_R).^((alpha)^2) * (2+(1-alpha)*u.*initial_R).^(1-alpha^2));
vpr = diff (v,initial_R);
funcr = @(initial_R) vpr;
kmax = 10;
initial_R(1) = 0;
initial_R(2) = 0.0070;
y(2) =  v(initial_R(2));
ypr(2) = funcr (initial_R(2));
tol = 0.000001
for k = 2:0.01:kmax
    initial_R(k+1) = initial_R(k) - y(k)/ypr(k);
    y(k+1) = feval (v,initial_R(k+1));
    c = initial_R(k+1) - initial_R(k);
    d = subs(c, initial_R(k+1), initial_R(k));
    if (abs(double(d))) < tol
        disp('The Numerical has ended')
    end
    ypr(k+1) = feval(funcr, initial_R(k+1));

end

I would like to know how else can I check the tolerance level because everytime I run this code with no matter how many variations, the if line throws an error.

Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for
variables.
Error in sym/double (line 868)
        Xstr = mupadmex('symobj::double', S.s, 0);
Related documentation

Answer 1

[Torsten]

syms initial_R
u = -100;
alpha = 2;
v = initial_R*exp((2*(1-alpha))/(2+(1-alpha)*u*initial_R)/(2-alpha*u*initial_R)^(alpha)^2*(2+(1-alpha)*u*initial_R)^(1-alpha^2))

v = 

vpr = diff(v,initial_R);
y = matlabFunction(v)

y = function_handle with value:
    @(initial_R)initial_R.*exp(1.0./(initial_R.*1.0e+2+2.0).^4.*1.0./(initial_R.*2.0e+2+2.0).^4.*-2.0)

ypr = matlabFunction(vpr)

ypr = function_handle with value:
    @(initial_R)exp(1.0./(initial_R.*1.0e+2+2.0).^4.*1.0./(initial_R.*2.0e+2+2.0).^4.*-2.0)+initial_R.*exp(1.0./(initial_R.*1.0e+2+2.0).^4.*1.0./(initial_R.*2.0e+2+2.0).^4.*-2.0).*(1.0./(initial_R.*1.0e+2+2.0).^4.*1.0./(initial_R.*2.0e+2+2.0).^5.*1.6e+3+1.0./(initial_R.*1.0e+2+2.0).^5.*1.0./(initial_R.*2.0e+2+2.0).^4.*8.0e+2)

kmax = 50;
initial_R_old = 1.0;
tol = 0.000001;
k = 0;
dx = 1.0;
while dx > tol && k < kmax
    k = k+1;
    initial_R_new = initial_R_old - y(initial_R_old)/ypr(initial_R_old);
    dx = abs(initial_R_old - initial_R_new);
    initial_R_old = initial_R_new;
end
initial_R_new

initial_R_new = 0

y(initial_R_new)

ans = 0