%Script file to show a sample solution curve
% to the geodesic on a torus.  The function
% torusDE.m needs to be written for the ODE
% solver.

%Solve the BVP:
[t,X]=ode23s(@torusDE,[0:0.01:10],[0,1,0,1/2]);

a=X(:,1);
b=X(:,3);

%This is a plot in the (alpha, beta) plane (beta
% is the "independent" variable)
figure(1)
plot(b,a);
%Alternatively, you could do two plots:
figure(2)
plot(t,a,t,b);

%To see the curve on the torus, we can use the following
%  set of commands:
figure(3)
[u,v]=meshgrid(linspace(0,2*pi,36),linspace(0,2*pi,36));
xx=(2+cos(u)).*cos(v);
yy=(2+cos(u)).*sin(v);
zz=sin(u);

x=(2+cos(a)).*cos(b);
y=(2+cos(a)).*sin(b);
z=sin(a);

A1=plot3(x,y,z,'k-');
set(A1,'LineWidth',3);
hold on
surf(xx,yy,zz);
shading interp
colormap(summer)
axis equal
hold off