function doall=glider
clear all; close all;
% Eli Tziperman: example for course in nonlinear dynamics and
% chaos:  Solve glider problem: v=velocity, u=angle:
% \dot{v}=-sin u -Dv^2; \dot{u}=(-cos u+v^2)/v
% also add equations for the horizontal and vertical coordinates (x,y):
% dx/dt=v*cos u; dy/dt=v*sin u;

% define parameters:
% ------------------
% friction coefficient:
D=0;
% time step: needs to be small enough for a given problem:
delta_t=0.01;
% Print (and save) every nprint steps:
nprint=1;
% Integration time:
T=15.0;

% Number of steps:
N=floor(T/delta_t);
NP=floor(N/nprint);

% Initial values:
u0=0; v0=0.4; x0=0; y0=0;
x(1:4,1)=[v0,u0,x0,y0]';

% Integrate:
for n=1:1:NP
    xp=x(:,n)';
    for j=1:1:nprint
      x_tilde=xp+delta_t*rhs(D,xp);
      xp=xp+0.5*( rhs(D,xp) + rhs(D,x_tilde) )*delta_t;
    end
    x(:,n+1)=xp';
    time(n)=n*delta_t;
end

% Plot simple time series:
plot(x(3,1:n),x(4,1:n))
hold on
xlabel('x')
ylabel('y')

% Quiver plot:
N=20;
Lu=pi*3;
Lv=3;
udot=zeros(N,N);
vdot=zeros(N,N);
[v,u] = meshgrid(0.2:Lv/N:Lv+0.2,0.2:Lu/N:Lu+0.2);

for i=1:N+1
    for j=1:N+1
        vdot(i,j)=-sin(u(i,j));
        udot(i,j)=(-cos(u(i,j))+v(i,j)^2)/v(i,j);
    end
end

figure
quiver(u,v,udot,vdot); xlabel('u'); ylabel('v'); title('phase plane')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ff=rhs(D,x)
% return the rhs of the dynamical system:

v=x(1);
u=x(2);

% dv/dt:
ff(1)=-sin(u) -D*v^2;
% du/dt:
ff(2)=(-cos(u)+v^2)/v;
% dx/dt:
ff(3)=v*cos(u);
% dy/dt:
ff(4)=v*sin(u);
