Method of Characteristics for Nozzle: Error with syms function

I am trying to create a Method of characteristics Code from JD Andersons book "Modern Compressible Flow with Historical Perspective" using example 11.1. Error arises when I try to obtain the coordinates of the points 2-7; where 7 is the n_div which I have taken.
Code:
clear;
clc;
close all;
%% Initial parameters %%
radian_to_deg = 180/pi;
degree_to_radian = pi/180;
prompt = "Enter a Mach Number ranging from 1.5 to 4: ";
Me = input(prompt);
prompt_1 = "Enter the number of divisions: ";
n_div = input(prompt_1);
%prompt_2 = "Enter the throat radius: ";
%TR = input(prompt_2);
g = 1.4;
mach_angle= [];
theta = [];
theta_a = [];
y=[];
x=[];
v=[];
v_a=[];
dydx_minus = [];
dydx_plus = [];
new_Mach=[];
%% get mach angle
ma_f = @(x) (asind(1/x)) ;
%% get theta_max
A = sqrt((g+1)/(g-1));
B = (g-1)/(g+1);
pm_f = @(x) (A*atand(sqrt(B*(x^2-1))) – atand(sqrt(x^2-1)));
theta_max = pm_f(Me)/2;
%% extraction of fractional part of theta_max
theta_i = sign(theta_max)*(abs(theta_max) – floor(abs(theta_max)));
%% Storing different values of theta_a
for i = 1:n_div
theta_a(i,1) = theta_i ;
theta_i = theta_i + 3.0;
end
v_a = theta_a;
%% Getting properties of a-1 characteristic line
ya=1;
xa=0;
%theta_a(1,1) = theta_i(1,1);
v_a(1,1) = theta_a(1,1);
%point 1
y(1,1) = 0;
Km_1 = theta_a(1,1) + v_a(1,1);
theta(1,1)= 0.5*(Km_1);
v(1,1) = 0.5*(Km_1);
%% get mach number for v(1,1)
syms x
f = @(x) (A*atand(sqrt(B*(x^2-1))) – atand(sqrt(x^2-1))-v(1,1));
f1 = matlabFunction(diff(f(x)));
%NR Method using a while loop
change = 10;
Mguess = 1.1;
while (change>1e-6)
Mnew = Mguess – f(Mguess)/f1(Mguess);
change = abs(Mguess – Mnew);
Mguess = Mnew;
end
new_Mach(1,1) = Mnew;
mach_angle(1,1) = ma_f(Mnew);
x(1,1) = -ya/tand(theta(1,1)-mach_angle(1,1));
%% rest of points:
for i = 2:n_div
Km = (theta_a(i,1)) + (v_a(i,1));
Kp = theta(i-1,1) – v(i-1,1);
theta(i,1) = 0.5*(Km+Kp);
v(i,1) = 0.5*(Km-Kp);

end
%% Calculating the Unknown Mach Number and corresponding Mach Angle using NR Method
syms x
for i=2:n_div
f = @(x) (A*atand(sqrt(B*(x^2-1))) – atand(sqrt(x^2-1))-v(i,1));
f1 = matlabFunction(diff(f(x)));
Mguess = 1.1;
change = 10;
while (change>1e-6)
Mnew = Mguess – f(Mguess)/f1(Mguess);
change = abs(Mguess – Mnew);
Mguess = Mnew;
end
new_Mach(i,1) = Mnew;
mach_angle(i,1) = ma_f(Mnew);
end

%% slopes and locations of points from 2-7
for i=2:n_div
dydx_minus(i,1) = tand(0.5*(theta(i)+theta(i)) – 0.5*(mach_angle(i)+mach_angle(i)));
dydx_plus(i,1) = tand(0.5*(theta(i) + theta(i-1))) + (0.5*(mach_angle(i) + mach_angle(i-1)));
x(i,1) = (y(i-1,1) – ya – (dydx_plus(i,1)*x(i-1)))/(dydx_minus(i,1) – dydx_plus(i,1));
y(i,1) = (dydx_plus(i,1)*(x(i,1)-x(i-1))) + y(i-1);
end
%slope and location of wall point 8
slopew1_minus = tand(theta_max);
slopew1_plus = dydx_plus(end,1);
xw = (y(end,1)-ya-(slopew1_plus*x(end,1)))/(slopew1_minus – slopew1_plus);
yw = (slopew1_plus*(xw(1,1)-x(end,1)))+y(end,1);
Error:
Error in MOC_prac (line 111)
y(i,1) = (dydx_plus(i,1)*(a_double(i,1)-a_double(i-1))) + y(i-1);

Unable to perform assignment because value of type ‘sym’ is not convertible to ‘double’.

Error in MOC_prac (line 111)
y(i,1) = (dydx_plus(i,1)*(x(i,1)-x(i-1))) + y(i-1);

Caused by:
Error using mupadengine/feval2char
Unable to convert expression containing symbolic variables into double array. Apply ‘subs’ function
first to substitute values for variables.
Is there any way I can resolve this error? Thank you!I am trying to create a Method of characteristics Code from JD Andersons book "Modern Compressible Flow with Historical Perspective" using example 11.1. Error arises when I try to obtain the coordinates of the points 2-7; where 7 is the n_div which I have taken.
Code:
clear;
clc;
close all;
%% Initial parameters %%
radian_to_deg = 180/pi;
degree_to_radian = pi/180;
prompt = "Enter a Mach Number ranging from 1.5 to 4: ";
Me = input(prompt);
prompt_1 = "Enter the number of divisions: ";
n_div = input(prompt_1);
%prompt_2 = "Enter the throat radius: ";
%TR = input(prompt_2);
g = 1.4;
mach_angle= [];
theta = [];
theta_a = [];
y=[];
x=[];
v=[];
v_a=[];
dydx_minus = [];
dydx_plus = [];
new_Mach=[];
%% get mach angle
ma_f = @(x) (asind(1/x)) ;
%% get theta_max
A = sqrt((g+1)/(g-1));
B = (g-1)/(g+1);
pm_f = @(x) (A*atand(sqrt(B*(x^2-1))) – atand(sqrt(x^2-1)));
theta_max = pm_f(Me)/2;
%% extraction of fractional part of theta_max
theta_i = sign(theta_max)*(abs(theta_max) – floor(abs(theta_max)));
%% Storing different values of theta_a
for i = 1:n_div
theta_a(i,1) = theta_i ;
theta_i = theta_i + 3.0;
end
v_a = theta_a;
%% Getting properties of a-1 characteristic line
ya=1;
xa=0;
%theta_a(1,1) = theta_i(1,1);
v_a(1,1) = theta_a(1,1);
%point 1
y(1,1) = 0;
Km_1 = theta_a(1,1) + v_a(1,1);
theta(1,1)= 0.5*(Km_1);
v(1,1) = 0.5*(Km_1);
%% get mach number for v(1,1)
syms x
f = @(x) (A*atand(sqrt(B*(x^2-1))) – atand(sqrt(x^2-1))-v(1,1));
f1 = matlabFunction(diff(f(x)));
%NR Method using a while loop
change = 10;
Mguess = 1.1;
while (change>1e-6)
Mnew = Mguess – f(Mguess)/f1(Mguess);
change = abs(Mguess – Mnew);
Mguess = Mnew;
end
new_Mach(1,1) = Mnew;
mach_angle(1,1) = ma_f(Mnew);
x(1,1) = -ya/tand(theta(1,1)-mach_angle(1,1));
%% rest of points:
for i = 2:n_div
Km = (theta_a(i,1)) + (v_a(i,1));
Kp = theta(i-1,1) – v(i-1,1);
theta(i,1) = 0.5*(Km+Kp);
v(i,1) = 0.5*(Km-Kp);

Unable to perform assignment because value of type ‘sym’ is not convertible to ‘double’.

Error in MOC_prac (line 111)
y(i,1) = (dydx_plus(i,1)*(x(i,1)-x(i-1))) + y(i-1);

Caused by:
Error using mupadengine/feval2char
Unable to convert expression containing symbolic variables into double array. Apply ‘subs’ function
first to substitute values for variables.
Is there any way I can resolve this error? Thank you! I am trying to create a Method of characteristics Code from JD Andersons book "Modern Compressible Flow with Historical Perspective" using example 11.1. Error arises when I try to obtain the coordinates of the points 2-7; where 7 is the n_div which I have taken.
Code:
clear;
clc;
close all;
%% Initial parameters %%
radian_to_deg = 180/pi;
degree_to_radian = pi/180;
prompt = "Enter a Mach Number ranging from 1.5 to 4: ";
Me = input(prompt);
prompt_1 = "Enter the number of divisions: ";
n_div = input(prompt_1);
%prompt_2 = "Enter the throat radius: ";
%TR = input(prompt_2);
g = 1.4;
mach_angle= [];
theta = [];
theta_a = [];
y=[];
x=[];
v=[];
v_a=[];
dydx_minus = [];
dydx_plus = [];
new_Mach=[];
%% get mach angle
ma_f = @(x) (asind(1/x)) ;
%% get theta_max
A = sqrt((g+1)/(g-1));
B = (g-1)/(g+1);
pm_f = @(x) (A*atand(sqrt(B*(x^2-1))) – atand(sqrt(x^2-1)));
theta_max = pm_f(Me)/2;
%% extraction of fractional part of theta_max
theta_i = sign(theta_max)*(abs(theta_max) – floor(abs(theta_max)));
%% Storing different values of theta_a
for i = 1:n_div
theta_a(i,1) = theta_i ;
theta_i = theta_i + 3.0;
end
v_a = theta_a;
%% Getting properties of a-1 characteristic line
ya=1;
xa=0;
%theta_a(1,1) = theta_i(1,1);
v_a(1,1) = theta_a(1,1);
%point 1
y(1,1) = 0;
Km_1 = theta_a(1,1) + v_a(1,1);
theta(1,1)= 0.5*(Km_1);
v(1,1) = 0.5*(Km_1);
%% get mach number for v(1,1)
syms x
f = @(x) (A*atand(sqrt(B*(x^2-1))) – atand(sqrt(x^2-1))-v(1,1));
f1 = matlabFunction(diff(f(x)));
%NR Method using a while loop
change = 10;
Mguess = 1.1;
while (change>1e-6)
Mnew = Mguess – f(Mguess)/f1(Mguess);
change = abs(Mguess – Mnew);
Mguess = Mnew;
end
new_Mach(1,1) = Mnew;
mach_angle(1,1) = ma_f(Mnew);
x(1,1) = -ya/tand(theta(1,1)-mach_angle(1,1));
%% rest of points:
for i = 2:n_div
Km = (theta_a(i,1)) + (v_a(i,1));
Kp = theta(i-1,1) – v(i-1,1);
theta(i,1) = 0.5*(Km+Kp);
v(i,1) = 0.5*(Km-Kp);

Unable to perform assignment because value of type ‘sym’ is not convertible to ‘double’.

Error in MOC_prac (line 111)
y(i,1) = (dydx_plus(i,1)*(x(i,1)-x(i-1))) + y(i-1);

Cart

Cart