a = 2;
b = 2;
T = 2*pi*sqrt(a^2+b^2);
kappa = abs(a)/(a^2+b^2);
nu = 1;
tau = b/(a^2+b^2);
kappa = (s)kappa;
nu = (s)nu;
tau = (s)tau;
A = (s)[ 0 kappa(s)*nu(s) 0
-kappa(s)*nu(s) 0 tau(s)*nu(s)
0 -tau(s)*nu(s) 0]
f = (s,y)reshape(A(s)*reshape(y,[3 3]),9,1); % Handle matrix ODE
sspan = linspace(0,T,20);
y0 = eye(3); % [e1 e2 e3;n1 n2 n3;b1 b2 b3]
[s,y] = ode45(f,sspan,y0(:));
y = reshape(y.,[3 3 numel(y)/9])
% Three basis vectors as functions of arc length
E = y(:,[1 4 7]).'; % [e11 ... e1m;e21 ... e2m;e31 ... e3m]
N = y(:,[2 5 8]).'; % [n11 ... n1m;n21 ... n2m;n31 ... n3m]
B = y(:,[3 6 9]).'; % [b11 ... b1m;b21 ... b2m;e31 ... b3m]
% Circular helix
t = s.'/sqrt(a^2+b^2);
x = a*cos(t);
y = a*sin(t);
z = b*t;
figure
plot3(x,y,z,'r')
hold on
axis equal
grid on