Trouble plotting solution to two 2nd order ODEs using forward Euler method

I have to graph the solution to two two-dimensional projectile motion including air resistance (drag) 2nd order ODEs that have been broken down to four 1st order ODEs using the forward Euler method and ode45.

The ODEs are:

This is what I have so far but I'm having trouble with the Euler solution:

clc; clear;

v0=85;%m/s
angle=51*pi/180;%rad
g=9.81;%m/s^2
m=0.145;%kg
r=0.036;%m
A=pi*r^2;%cross section area
x0=0;y0=0.9;%initial position
%Case 1
C=0.098;
rho=1.225;%kg/m^3
[x1,y1,v1,v2]=projectile_motion(v0,angle,g,m,A,C,rho,y0,x0);%Euler solution
%ode45 solution
ts=0:0.1:11;
ic=[x0,y0,v0*cos(angle),v0*sin(angle)];%initial conditions
[t,x] = ode45(@f,ts,ic);
%Analytical without drag
tse=0:0.1:13.45;%time span extended to see full trajectory for analy. part
xa=x0+v0*cos(angle)*tse;
ya=y0+v0*sin(angle)*tse-0.5*g*tse.^2;



figure(1)%plot of x v/s y
hold on
plot(x1,y1,'*r',x(:,1),x(:,2),'ob');
plot(xa,ya,'.k',x1,y1,'r');
plot(x(:,1),x(:,2),'b',xa,ya,'k')
hold off
xlabel('x(m)'); ylabel('y(m)');
title('2D path of the projectile');
legend('Forward Euler','ode45','without drag');
ylim([-10,240]);

figure(2)
subplot(2,1,1)
hold on
plot(ts,v1(:,1),'*r',t,x(:,3),'ob');
plot(ts,v1(:,1),'r',t,x(:,3),'b');
hold off
xlabel('t (s)'); ylabel('v_x (m/s)');
title('v_x vs time'); legend('Forward Euler','ode45');
subplot(2,1,2)
hold on
plot(ts,v2(:,1),'*r',t,x(:,4),'ob');
plot(ts,v2(:,1),'r',t,x(:,4),'b');
hold off
xlabel('time (seconds)')
ylabel('v_y (m/s)')
legend('Euler','ode45')
title('v_y vs time')
%local function for ode45

%forward euler method
for n=1:1:length(t)-1
    v(n)=(vx(n)^2+vy(n)^2)^0.5;
    ay(n+1)=-g-(D/m)*v(n)*vy(n);
    vy(n+1)=vy(n)+ay(n)*dt;
    y(n+1)=y(n)+vy(n)*dt+0.5*ay(n)*dt^2;
    ax(n+1)=-(D/m)*v(n)*vx(n);
    vx(n+1)=vx(n)+ax(n)*dt;
    x(n+1)=x(n)+vx(n)*dt+0.5*ax(n)*dt^2;
    
end

Function "projectile_motion":
function [x,y,vx,vy]=projectile_motion(v0,angle,g,m,A,C,rho,y0,x0)% Solution by Euler method
D=0.5*(rho*C*A);
%initail conditions
vx0=v0*cos(angle);
vy0=v0*sin(angle);
ax0=-(D/m)*v0*vx0;
ay0=-g-(D/m)*v0*vy0;
dt=0.1;%step size
end_t=11;%final time
t=0:dt:end_t;%time span
%allocating arrays 
x=zeros(1,length(t));
y=zeros(1,length(t));
vx=zeros(length(t));
vy=zeros(length(t));
v=zeros(length(t));
ax=zeros(length(t));
ay=zeros(length(t));
x(1)=x0;
y(1)=y0;
vx(1)=vx0;
vy(1)=vy0;
ax(1)=ax0;
ay(1)=ay0;
end


function dtdx=f(t,x)%system of DEs
    g=9.81; m=0.145; r=0.036; A=pi*r^2;
    rho=1.225; C=0.098; D=0.5*(rho*C*A);
    v=(x(3)^2+x(4)^2)^0.5;%magnitude of velocity.
    dtdx=[x(3);x(4);-(D/m)*v*x(3);-g-(D/m)*v*x(4)];
end

Answer 1

[Alan Stevens]

The following sorts out your Euler integration:

v0=85;%m/s
angle=51*pi/180;%rad
g=9.81;%m/s^2
m=0.145;%kg
r=0.036;%m
A=pi*r^2;%cross section area
x0=0;y0=0.9;%initial position
%Case 1
C=0.098;
rho=1.225;%kg/m^3
[x1,y1,v1,v2]=projectile_motion(v0,angle,g,m,A,C,rho,y0,x0);%Euler solution
%ode45 solution
ts=0:0.1:11;
ic=[x0,y0,v0*cos(angle),v0*sin(angle)];%initial conditions
[t,x] = ode45(@f,ts,ic);
%Analytical without drag
tse=0:0.1:13.45;%time span extended to see full trajectory for analy. part
xa=x0+v0*cos(angle)*tse;
ya=y0+v0*sin(angle)*tse-0.5*g*tse.^2;



figure(1)%plot of x v/s y
hold on
plot(x1,y1,'*r',x(:,1),x(:,2),'ob',xa,ya,'.k');
hold off
xlabel('x(m)'); ylabel('y(m)');
title('2D path of the projectile');
legend('Forward Euler','ode45','without drag');
ylim([-10,240]);

figure(2)
subplot(2,1,1)
hold on
plot(ts,v1,'*r',t,x(:,3),'ob');
% plot(ts,v1,'r',t,x(:,3),'b');
hold off
xlabel('t (s)'); ylabel('v_x (m/s)');
title('v_x vs time'); legend('Forward Euler','ode45');
subplot(2,1,2)
hold on
plot(ts,v2,'*r',t,x(:,4),'ob');
hold off
xlabel('time (seconds)')
ylabel('v_y (m/s)')
legend('Euler','ode45')
title('v_y vs time')


% Function "projectile_motion":
function [x,y,vx,vy]=projectile_motion(v0,angle,g,m,A,C,rho,y0,x0)% Solution by Euler method
D=0.5*(rho*C*A);
%initial conditions
vx0=v0*cos(angle);
vy0=v0*sin(angle);
dt=0.1;%step size
end_t=11;%final time
t=0:dt:end_t;%time span
%allocating arrays 
xy=zeros(2,length(t));
vxvy=zeros(2,length(t));
xy(:,1)=[x0; y0];
vxvy(:,1)=[vx0; vy0];

% Simple Euler updating
for i=1:numel(t)-1
  
    xy(:,i+1) = xy(:,i) + vxvy(:,i)*dt;  % update position components
    vx = vxvy(1,i); vy = vxvy(2,i);
    v = sqrt(vx^2 + vy^2);  % current velocity
    dvdt = [-(D/m)*v*vx;-g-(D/m)*v*vy]; % current acceleration components
    vxvy(:,i+1) = vxvy(:,i) + dvdt*dt;   % update velocities
end

% Extract individual component positions and velocities
x = xy(1,:);
y = xy(2,:);
vx = vxvy(1,:);
vy = vxvy(2,:);

end


function dtdx=f(~,x)%system of DEs
    g=9.81; m=0.145; r=0.036; A=pi*r^2;
    rho=1.225; C=0.098; D=0.5*(rho*C*A);
    v=(x(3)^2+x(4)^2)^0.5;%magnitude of velocity.
    dtdx=[x(3);x(4);-(D/m)*v*x(3);-g-(D/m)*v*x(4)];
end