Projectile 2D plot with drag using ODE45

I am doing this interesting project to plot a 2D Trajectory of projectile under an air drag. I have written two functions for that, function f.m where I have listed the ODEs describing the motion and function QuadDrag.m where I use ODE45 fuction of matlab and plot the graph.

Both of them look fine, but I am getting an error again and again. Moreover I want to stop the plot when y value becomes 0 using while loop. But I cant seem to get the correct plot and errors.

Here is my code below

function [xp] = f(t,z)
%====================================================================
% #Program: f.m
% #Purpose: define equations of motion for projectile motion
%           with quadratic drag in 'state variable form'.
%====================================================================
 
global m c C A 
m=0.005;        % mass of the projectile
g=9.81;         % acceleration due to gravity constant

t(1)=0;         %initial time

rho = 1.173;    %Density in kg/m^3
c=(rho*A*C)/2;  %assuming small letter 'c' to be the combined constant

xp=zeros(4,1);  %matrix xp

%-- Equations of Motion (State Variable Form) --
xp(1)=z(2);                                          % x(velocity)
xp(2)=(-c/m)*sqrt( ((z(2))^2) + ((z(4))^2) )*z(2);   % x(acceleration)
xp(3)=z(4);                                          % y(velocity)
xp(4)=-g-(c/m)*sqrt( ((z(2))^2) + ((z(4))^2) )*z(4); % y(acceleration)
 
end

and QuadDrag.m

%QuadDRag.m
%-- Initial Conditions --
theta = 37.6; %in degrees
v_0 = 75;     %initial velocity in m/s

z1o=-2.1;              % x-(initial x position)
z2o= v_0*cosd(theta);  % x-(initial x velocity)
z3o=0;                 % y-(initial y position)
z4o= v_0*sind(theta);  % y-(initial y velocity)

y = [z1o;z2o;z3o;z4o];

toggle = true; 

while toggle
    
[t,z]=ode45('f',t,y);   % Solve DE With Quadratic Drag

 if y(3)<= 0            %Stop when projectile reaches the ground 
    toggle = false;     %i.e x axis when y value is equal to 0 
 end
end

%-- Plots --
plot(z(:,1),z(:,3));
title('x vs. y', 'FontSize', 14)
xlabel('x', 'FontSize', 14);
ylabel('y', 'FontSize', 14);
legend('Quadratic Drag');

Please tell me if you find any mistakes or suggest changes to improve the same.

Thank you so much for your time.

[madhan ravi]

can you upload the latex form of the equation?

[Prathamesh Toraskar]

x'' = -\frac{c}{m}*\sqrt{(x'^{2}+y'^{2})} *x'

y'' = -\frac{c}{m}*\sqrt{(x'^{2}+y'^{2})} *y'

[Prathamesh Toraskar]

Those are the equations! where x is x position and y is y position. Hence x' = x Velocity and y' = y Velocity and x'' = x accelerationaand y'' = y acceleration

[madhan ravi]

I suspect you should be using bvpc4 also respond to Steven Lord's answer

Answer 1

[Steven Lord]

Take a look at the ballode example. It uses an events function to stop the ODE solver when the bouncing ball hits the floor. You may be able to adapt it to your system of ODEs.

A few other suggestions:

1. Don't specify the function as the first input to ode45 using a char vector. Use a function handle (which could be an anonymous function) instead.
2. Don't use global variables. We can't run your code because we don't know what values those global variables have, and anyone with access to the global workspace could affect your code's computation without you knowing. [How do you know that sqrt or some function called by ode45 doesn't set a global variable C or A?] See the ode45 documentation page for instructions on how to parameterize your ODE function.
3. If you want to see the position of the projectile while the ODE solver is running, consider using an OutputFcn. See the documentation for the odeset function for a description of this option.
4. When you post a question where you're receiving an error that you want to eliminate, please post the full text (everything displayed in red) of the error message you receive. That usually contains useful information that will help others determine the cause of the problem more quickly.

[Prathamesh Toraskar]

Thank you for the suggestions. I update the question and respond ASAP.

Answer 2

[madhan ravi]

https://www.mathworks.com/matlabcentral/answers/15119-use-of-ode45-for-projectile-trajectory-including-drag#answer_20524