maximaで遊んだ時のコード

animation.mac

with_slider_draw(k,makelist(i,i,1,7),
	proportional_axes = xy,
	xrange = [-3,3],
	yrange = [-3,3],
	xaxis=true,
	yaxis=true,
	grid = true,
	title = "animation" ,
	label([sconcat("k=",k-4),-1,2]),
	point_type=7,
	points([[k-4, k-4]])
);

curve.mac

f(p, p0, v0, g, t0, t1, dt) := block([sol, s, t, s1, p1: diff(p(s), s), p2: diff(p(s), s, 2)], sol: rk([s1, ((s1^2)*(p1.p2)+g*p1[3])/(p1.p1)], [s, s1], [p0, v0], [t, t0, t1, dt]), makelist([l[1], l[2], p(l[2]), vector(p(l[2]), l[3]^2*subst(l[2], s, p2-p2.p1/(p1.p1)*p1)), vector(p(l[2]), l[3]*subst(l[2], s, p1))], l, sol))$

load("draw");

/*
sol: f(lambda([s], [s, 0, s^2]), -1, 0, -10, 0, 5, 0.1)$
draw3d (proportional_axes=xyz, xrange=[-5, 5], yrange=[-5, 5], head_length=0.1,makelist(p[4], p, sol));
*/

sol: f(lambda([s], [cos(s), sin(s), s]), 0, 0, -10, 0, 1, 0.1)$
draw3d (proportional_axes=xyz,head_length=0.5,makelist(p[4], p, sol), parametric(cos(s), sin(s), s, s, -2*%pi, 0));

surface.mac

f(p, iuv, g, t0, t1, dt) := block(
	[sol, u, v, t, u1, v1, pu: diff(p(u, v), u), pv: diff(p(u, v), v), puu: diff(p(u, v), u, 2), puv: diff(p(u, v), u, 1, v, 1), pvv: diff(p(u, v), v, 2), q, m, r],
	q: u1^2*puu+2*u1*v1*puv+v1^2*pvv,
	m: matrix([pv.pv, -pu.pv], [-pu.pv, pu.pu])/((pu.pu)*(pv.pv)-(pu.pv)^2),
	r: -m.([pu.q, pv.q]+g*[pu[3], pv[3]]),
	sol: rk([u1, v1, r[1][1], r[2][1]], [u, v, u1, v1], iuv, [t, t0, t1, dt]),
	makelist([l[1], l[2], l[3], p(l[2], l[3])], l, sol)
)$

/*
sol: f(lambda([u, v], [u, v, u^2+v^2]), [1, 1, 0, 0], 10, 0, 3, 0.1)$

draw3d(proportional_axes=xyz, point_size=5, point_type=7, points(makelist(p[4], p, sol)), parametric_surface(u, v, u^2+v^2, u, -1, 1, v, -1, 1));
*/

sol: f(lambda([r, t], [r*cos(t), r*sin(t), r+t]), [1, 0, 0, 0], 10, 0, 5, 0.1)$
draw3d(proportional_axes=xyz, point_size=3, point_type=7, points(makelist(p[4], p, sol)), color=green, parametric_surface(r*cos(t), r*sin(t), r+t, r, 0, 5, t, -5*%pi, 0));

with_slider_draw3d(
	k,makelist(k,k,1,50),
	proportional_axes=xyz,
	dimensions=[500,500],
	xrange=[-10,10],
	yrange=[-10,10],
	zrange=[-15, 5],
	point_type=7,
	points([sol[k][4]]),
	color=green,
	parametric_surface(r*cos(t), r*sin(t), r+t, r, 0, 10, t, -5*%pi, 0)
)$