clear
epsilon=7.5e-4;
delta=.0124;
g=9.8;

a=0;
b=.1;
N=1000;
h=(b-a)/N;

%f = inline('-g+delta*(v/r).^(1.5) +epsilon*v.^2/r','v','r');
%g = inline('epsilon/3*r','v','r');
f=@(v,r) g-delta*(v/r).^(1.5) -epsilon*v^2/r;
g=@(v,r) epsilon/3*r;

v(1)=0;
r(1)=.0001;
t(1)=0;

for j = 2 : N+1
    k1 = h * f(v(j-1),r(j-1));
    m1 = h * g(v(j-1),r(j-1));
    k2 = h * f( v(j-1) + k1/2, r(j-1) + m1/2 );
    m2 = h * g( v(j-1) + k1/2, r(j-1) + m1/2 );
    k3 = h * f( v(j-1) + k2/2, r(j-1) + m2/2 );
    m3 = h * g( v(j-1) + k2/2, r(j-1) + m2/2 );
    k4 = h * f( v(j-1) + k3, r(j-1) + m3 );
    m4 = h * g( v(j-1) + k3, r(j-1) + m3 );
    v(j) = v(j-1) + (k1 + 2*k2 + 2*k3 + k4) / 6;
    r(j) = r(j-1) + (m1 + 2*m2 + 2*m3 + m4) / 6;
    t(j)=t(j-1)+h;
  end