Hi, I'm trying to integrate the Taylor-Maccoll equation in order to determine the velocity flowfield in a conical shockwave, moving at hypersonic speed. As explained in the photo ive attached, the equation can be represented as a standard ODE with solution vector. I have initial conditions for both vr and vr', which I am then supposed to place into the second vector, y'. I then integrate y' using the Runge-Kutta method and at some point in the integration vr' should equal (roughly) 0. I have a step size of h, which defines the increment in theta (polar coordinates) I take from the shockwave, towards the cone angle.
I am aware of two codes which have been uploaded to the MATLAB servers which I can download, however these codes do not include the discrete values of velocity at every angle of theta, and so I must code these a different way.
In my code I have manually typed up the Runge-Kutta method but I think there is some problem in my code as my vrdash values are not equating to zero at the correct theta value (noweher near in fact!).
vr = 0.9078;
vrdash = -0.1262;
theta_s = 14.35;
y1 = zeros(1, length(x));
y2 = zeros(1, length(x));
y1(1) = vr;
y2(1) = vrdash;
h = -0.01;
x = theta_s:h:0;
g = 1.4;
t(1) = 0;
A = (g-1)/2;
F1 = @(t, y1, y2) y2;
F2 = @(t, y1, y2) (y1*y2^2-A*(1-y1^2-y2^2)*(2*y1+y2*cotd(theta_s)))...
/(A*(1-y1^2-y2^2) - y2^2);
k1 = h*F1(t, y1(i), y2(i));
m1 = h*F2(t, y1(i), y2(i));
k2 = h*F1(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*k1);
m2 = h*F2(t+0.5*h, y1(i)+0.5*h, y2(i)+0.5*h*m1);
k3 = h*F1(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*k2));
m3 = h*F2(t+0.5*h, (y1(i)+0.5*h), (y2(i)+0.5*h*m2));
k4 = h*F1(t+h, (y1(i)+h),(y2(i)+k3*h));
m4 = h*F2(t+h, (y1(i)+h),(y2(i)+m3*h));
y1(i+1) = y1(i) + (1/6)*(k1+2*k2+2*k3+k4)*h;
y2(i+1) = y2(i) + (1/6)*(m1+2*m2+2*m3+m4)*h;
I'm sorry if I havent explained this properly, if there's anything i'm missing I'd be happy to explain. Any help would be greatly appreciated as I've been stuck on this code for nearly a month now. Thanks!