function prel2q3()   
t_end = 100; x0=0; y0=0; vx0=300; vy0 = 400; 
  m = 2;   c = 0.1;    g = 10;
p.m = m; p.c = c  ;  p.g = g;
tspan = [0 t_end];  n = 20000;
z0    = [x0 y0 vx0 vy0]';


[tarray zarray] = rk2(tspan,z0,n,p);
% The next line nicely prints out zarray(end,1) 
disp(['x(100s)='  num2str(zarray(end,1))  'm'])
end


function [tmat zmat] = rk2(tspan,z0,n,p)
%Midpoint integration. No problems here. 
tmat   = linspace(tspan(1),tspan(2),n+1); 
h      = tmat(2)-tmat(1);
zmat   = zeros(n+1,length(z0)); 
zmat(1,:) = z0';
 for i=1:n;
    z           = zmat(i,:)';     t   = tmat(i);
    ztemp      =  z + (h/2)*rhs(t,z,p); 
    znew        = z + h*rhs((t+h/2),ztemp,p); 
    zmat(i+1,:) = znew';
 end
end


function  zdot = rhs(t,z,p)
m=p.m; c = p.c;  g = p.g;
r= z(1:2);  vvec = z(3:4); %position and velocity vectors
v = norm(vvec);  uv = vvec/v; % mag of vel, unit vector in dir. of vel.

Fdrag = -c*v^2*uv; % drag force, a 2-comp vector

rdot    = vvec;                   % first two ODEs
vvecdot = (1/m)*(Fdrag -[0 m*g]'); % 2nd   two ODEs (F = ma)

zdot = [rdot; vvecdot];  % All 4 comps of rhs.
end