Question: Making specific calculation computable

Hello,

I want to build a paper model of an immersion of real projective plane made out of paper strips. I want to use maple to plot the curves which represent the paper strips. I have written a code shown below which computes the curve for one of the strips. In this code, alpha is set to Pi/2 and the computation goes through. However the computation should be performed for various alpha. But if I set alpha=0, then f (see code below) is not computable, if I set alpha=Pi/4, the computation for DkNorm already takes ages. So here's the Q: How should I write the code, so that I get a plot for arbitrary alpha? Perhaps all computations numerically? But how?

restart; assume(x::real,y::real,t::real);
with(plots):with(LinearAlgebra):
M:=Re(1/(z^3-z^(-3)+sqrt(5))*<I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))>)+<0,0,2/5>:
B:=M/VectorNorm(M,2)^2:
BReal:=subs(z=x+I*y,B): BReal:=simplify(BReal) assuming real:
#plot3d([BReal(1),BReal(2),BReal(3)],x=-1..1,y=-1..1,grid=[80, 80]);
DxB:=<diff(BReal(1),x),diff(BReal(2),x),diff(BReal(3),x)>: DxB:=simplify(DxB) assuming real:
DyB:=<diff(BReal(1),y),diff(BReal(2),y),diff(BReal(3),y)>: DyB:=simplify(DyB) assuming real:
DxB0:=subs(x=0,y=0,DxB):
DyB0:=subs(x=0,y=0,DyB):


alpha:=Pi/2:
dir:=exp(I*alpha):
c:=dir*t:
cReal:=<Re(c),Im(c)>:
k:=subs(x=cReal(1),y=cReal(2),BReal): k:=simplify(k) assuming real:
#spacecurve([k(1),k(2),k(3)],t=0..1,color=black,axes=normal);
Dk:=<diff(k(1),t),diff(k(2),t),diff(k(3),t)>: Dk:=simplify(Dk) assuming real:
DkNorm:=VectorNorm(Dk,2): DkNorm:=simplify(DkNorm) assuming real:
DkN:=Dk/DkNorm: DkN:=simplify(DkN) assuming real:
DDkN:=<diff(DkN(1),t),diff(DkN(2),t),diff(DkN(3),t)>: DDkN:=simplify(DDkN) assuming real:
DxBc:=subs(x=cReal(1),y=cReal(2),DxB): DxBc:=simplify(DxBc) assuming real:
DyBc:=subs(x=cReal(1),y=cReal(2),DyB): DyBc:=simplify(DyBc) assuming real:
N:=CrossProduct(DxBc,DyBc): N:=simplify(N) assuming real:
NN:=N/VectorNorm(N,2): NN:=simplify(NN) assuming real:
q:=Matrix([[1+NN(1)*I,NN(2)+NN(3)*I],[-NN(2)+NN(3)*I,1-NN(1)*I]]):
x:=Matrix([[DkN(1)*I,DkN(2)+DkN(3)*I],[-DkN(2)+DkN(3)*I,-DkN(1)*I]]):
vv:=MatrixMatrixMultiply(q,MatrixMatrixMultiply(x,MatrixInverse(q))):
v:=<Im(vv[1,1]),Re(vv[1,2]),Im(vv[1,2])>: v:=simplify(v) assuming real:
f:=v(1)*DDkN(1)+v(2)*DDkN(2)+v(3)*DDkN(3): f:=simplify(f) assuming real:


dirV:=<Re(dir),Im(dir)>:
DBdirV:=<DxB0(1)*dirV(1)+DyB0(1)*dirV(2),DxB0(2)*dirV(1)+DyB0(2)*dirV(2),DxB0(3)*dirV(1)+DyB0(3)*dirV(2)>:
DBdirVN:=DBdirV/VectorNorm(DBdirV,2):
initial:=arccos(DBdirVN(1)):


fS:=subs(t=s,f):
DkNormS:=subs(t=s,DkNorm):
sys:=diff(theta(s),s)=fS, diff(b1(s),s)=cos(theta(s))*DkNormS, diff(b2(s),s)=sin(theta(s))*DkNormS:
p:=dsolve({sys,theta(0)=initial,b1(0)=0,b2(0)=0},{theta(s),b1(s),b2(s)},type=numeric):
plo:=odeplot(p,[b1(s),b2(s)],0..1,numpoints=200):
#display(plo);
t0:=0: t1:=0.2: t2:=0.4: t3:=0.6: t4:=0.8: t5:=1:
bx0:=subs(p(t0),b1(s)): bx1:=subs(p(t1),b1(s)): bx2:=subs(p(t2),b1(s)): bx3:=subs(p(t3),b1(s)): bx4:=subs(p(t4),b1(s)): bx5:=subs(p(t5),b1(s)):
by0:=subs(p(t0),b2(s)): by1:=subs(p(t1),b2(s)): by2:=subs(p(t2),b2(s)): by3:=subs(p(t3),b2(s)): by4:=subs(p(t4),b2(s)): by5:=subs(p(t5),b2(s)):
pts:=plot(Vector([bx0,bx1,bx2,bx3,bx4,bx5]), Vector([by0,by1,by2,by3,by4,by5]), style=point, symbol=solidcircle, color=blue):
display(plo,pts);