MATLAB: System of arc-length defined ODEs with ode45

odeode systemsode45

I am working on a code that has a system of ODEs, but I have never worked with systems with ode45. In the part of the code that I have included, S is the arclength (which is basically the time step of this problem), th is the angle (theta) of the graph, R is the x coordinate, and Z is the y coordinate.
When I run the program as shown below, I get simply a matrix full of NaN, even when I change the initial R value to 0.0001 or something.
Any help would be appreciated.
Also, what is the output? I only want to graph R and Z, not theta
function yp=program(S,y)
th=y(1);
R=y(2);
Z=y(3);
dthdS=-sin(th)/R+Z-2*H;
dRdS=cos(th);
dZdS=sin(th);
yp=[dthdS; dRdS; dZdS];
end
[S,Y]=ode45(@program, [0, 1], [0, 0, 0])

Best Answer

As with your 'How can I graph a "system" of ODEs?' post, I suggest:
        % dThXYdS(1) = Theta, dThXYdS(2) = R(t) = x,  dThXYdS(3) = Z(t) 
        H = 10;
        dThXYdS = @(t,ThXY) [-sin(ThXY(1))/ThXY(2) + ThXY(3) - 2*H;  cos(ThXY(1));  sin(ThXY(1))];
        x0 = [0.1; 0.1; 0];
        Tspan = [0:0.01:2]';
        [T ThXY] = ode45(dThXYdS, Tspan, x0);
        figure(8)
        plot(ThXY(:,2), ThXY(:,3))
        xlabel('R(S)')
        ylabel('Z(S)')
        grid
Except for the axis labels in the plot, I used the variable designation from your previous post (and my previous answer) rather than change it to match your current variable designation.
The reason you are getting a matrix of NaNs is that your initial conditions are [0 0 0]. So R is zero when theta is zero and of course sin(theta) will be zero as well. By convention, (0/0) = NaN. An initial NaN in a recursive calculation such as yours creates a matrix of NaNs.
Related Question