function main
%Set initial values for r, theta, z
r00 = 10;
theta00 = 0;
z00 = 1.0;
y0 = [r00 theta00 z00];
% Set model parameters
z0 = 0.0;
a = 0.2;
gamma_inf = 1000.0;
nu = 1.3;
% Set integration period
tspan = 0:0.001:10;
% Call integrator
fun = @(t,y)[-a*y(1);gamma_inf/(2*pi*y(1)^2)*(1-exp(-a*y(1)^2/(2*nu)));2*a*(y(3)-z0)];
[T Y] = ode45(fun,tspan,y0);
R = Y(:,1);
THETA = Y(:,2);
Z = -Y(:,3);
[x y] = pol2cart(THETA,R);
h=plot3(x,y,Z,'r', 'LineWidth', 1')
hold on;
plot3 ([0,0],[0,0],[Z(end),0],'--b')
grid on
xlabel x-axis
ylabel y-axis
zlabel z-axis
zticklabels({'0','10','20','30','40','50','60'})
set(h,'Zdir','reverse')
hold off
