summaryrefslogtreecommitdiff
path: root/texmf-dist/scripts/ketpic/ketlib/maximaL
diff options
context:
space:
mode:
Diffstat (limited to 'texmf-dist/scripts/ketpic/ketlib/maximaL')
-rw-r--r--texmf-dist/scripts/ketpic/ketlib/maximaL/fourier_sec.max260
-rw-r--r--texmf-dist/scripts/ketpic/ketlib/maximaL/matoperation.max47
-rw-r--r--texmf-dist/scripts/ketpic/ketlib/maximaL/poincare.mac71
-rw-r--r--texmf-dist/scripts/ketpic/ketlib/maximaL/rkfun.lisp99
4 files changed, 477 insertions, 0 deletions
diff --git a/texmf-dist/scripts/ketpic/ketlib/maximaL/fourier_sec.max b/texmf-dist/scripts/ketpic/ketlib/maximaL/fourier_sec.max
new file mode 100644
index 00000000..c70f7f46
--- /dev/null
+++ b/texmf-dist/scripts/ketpic/ketlib/maximaL/fourier_sec.max
@@ -0,0 +1,260 @@
+/*
+Package fourier_sec, to study piecewise defined functions
+(sectionally defined, hence the name).
+
+Author: Jose A. Vallejo
+ Faculty of Sciences
+ Universidad Autonoma de San Luis Potosi (Mexico)
+ http://galia.fc.uaslp.mx/~jvallejo
+
+It is assumed that the functions are defined using the format
+ if ... then ... else ...
+as, for example,
+(%i1) absolute(x):=if (x<=0) then -x elseif (x>0) then x$
+(it will not work with 'abs').
+It provides three functions:
+1. 'paritycheck', to check if a given piecewise defined
+ function is even or odd (or none). For instance,
+ (%i2) paritycheck(absolute(x),x);
+ (%o2) 0
+ A '0' means 'even', '1' is 'odd', and '-1' is 'none'
+2. 'fourier_sec_coeff', to compute the Fourier sine and cosine
+ coefficients. The function must have compact support.For
+ instance, for the unit step we have:
+ (%i3) h(x):=if (-%pi<=x and x<0) then 0 elseif (0<=x and x<=%pi) then 1$
+ (%i4) fourier_sec_coeff(h(x),x);
+ (%o4) [1/2,0,-((-1)^n-1)/(%pi*n)]
+3. 'fourier_sec_series', to compute the Fourier series, truncated
+ or not. The series for the unit step to the 10th order:
+ (%i5) fourier_sec_series(h(x),x,10);
+ (%o5) (2*sin(9*x))/(9*%pi)+(2*sin(7*x))/(7*%pi)+(2*sin(5*x))/(5*%pi)+(2*sin(3*x))/(3*%pi)+(2*sin(x))/%pi+1/2
+ and the whole Fourier series
+ (%i6) fourier_sec_series(h(x),x,inf);
+ (%o6) (sum(((1/n-(-1)^n/n)*sin(n*x))/%pi,n,1,inf))+1/2
+*/
+
+
+load(fourie)$
+
+load(simplify_sum)$
+
+paritycheck(funvar,x):=block([subintervals,subvalues,tmp1,tmp2,token1,token0,tmp3,
+ subvalues_list,subintervals_list,icentral:0,middle,minusmiddle,token2,
+ side_subintervals_list,L,side_subvalues_list,
+ zero_subintervals,tmp4,tmp5,non_zero_side_subvalues_list,
+ non_zero_side_subintervals_list,LL,token3,expr1,expr2],
+ local(M,N,P,Q,count1,count2,count3,count4),
+ subintervals:makelist(part(funvar,i),i,makelist(2*k-1,k,1,(length(funvar)-2)/2)),
+ subvalues:makelist(part(funvar,i),i,makelist(2*k,k,1,(length(funvar)-2)/2)),
+ for j:1 thru length(subintervals) do (if operatorp(subintervals[j],["<",">","<=",">="]) then tmp1[j]:1 else tmp1[j]:0),
+ tmp1:makelist(tmp1[j],j,1,length(subintervals)),
+ tmp2:sublist_indices(tmp1,lambda([x],x=1)),
+
+ /* if length(tmp2)=0 all the subintervals in the domain of funvar are bounded! */
+
+ if is(equal(length(tmp2),0)) then
+ (
+ subvalues_list:copylist(subvalues),
+ tmp3:copylist(subintervals),
+ for j:1 thru length(tmp3) do
+ (
+ M[j]:makelist(part(tmp3[j],k),k,1,length(tmp3[j])),
+ for m:1 thru 2 do N[j,m]:makelist(part(M[j],m,n),n,1,length(M[j])),
+ P[j]:append(N[j,1],N[j,2]),Q[j]:sort(delete(x,P[j]),"<")
+ ),
+ subintervals_list:makelist(Q[j],j,1,length(tmp3)),
+ for i:1 thru length(subintervals_list) do (if is(lmin(subintervals_list[i])*lmax(subintervals_list[i])<0) then icentral:i),
+ if is(icentral>0) then middle:subvalues_list[icentral],
+ if is(icentral>0) then
+ (if evenfunp(middle,x) then token2:0 elseif oddfunp(middle,x) then token2:1 else return(-1))
+ else
+
+ /* now we analyze what happens if icentral=0, so there are only side intervals */
+
+ (
+ /* as before, there must be an even number of symmetric intervals, otherwise -1 */
+ side_subintervals_list:copylist(subintervals_list),
+ if not(evenp(length(side_subintervals_list))) then return(-1),
+ L:length(side_subintervals_list)/2,
+ for k:1 thru L do count1[k]:charfun(is(equal(side_subintervals_list[k],reverse(map("-",side_subintervals_list[2*L+1-k]))))),
+ if sum(count1[j],j,1,L)#L then return(-1),
+ side_subvalues_list:copylist(subvalues_list),
+ for k:1 thru L do count2[k]:charfun(is(equalp(side_subvalues_list[k],ratsubst(-x,x,side_subvalues_list[2*L+1-k])))),
+ if is(equal(sum(count2[j],j,1,L),L)) then token3:0,
+ for k:1 thru L do count3[k]:charfun(is(equalp(side_subvalues_list[k],-ratsubst(-x,x,side_subvalues_list[2*L+1-k])))),
+ if is(equal(sum(count3[j],j,1,L),L)) then token3:1,
+ if is(not(equalp(token3,0)) and not(equalp(token3,1))) then return(-1) ,
+ return(token3)
+ ),
+
+ /* and what happens if icentral#0 */
+
+ side_subintervals_list:delete(subintervals_list[icentral],subintervals_list),
+ if is(equal(length(side_subintervals_list),0)) then (if is(equal(token1,token2)) then return(token1*token2) else return(-1)),
+ if not(evenp(length(side_subintervals_list))) then return(-1),
+ L:length(side_subintervals_list)/2,
+ for k:1 thru L do count1[k]:charfun(is(equal(side_subintervals_list[k],reverse(map("-",side_subintervals_list[2*L+1-k]))))),
+ if sum(count1[j],j,1,L)#L then return(-1),
+ side_subvalues_list:delete(subvalues_list[icentral],subvalues_list),
+
+ /* remove those subintervals in which funvar vanishes */
+
+ zero_subintervals:sublist_indices(side_subvalues_list,lambda([x],x=0)),
+ tmp4:copylist(side_subintervals_list),
+ for j:1 thru length(zero_subintervals) do
+ (tmp4:delete(side_subintervals_list[zero_subintervals[j]],tmp4)),
+ non_zero_side_subintervals_list:tmp4,
+ tmp5:copylist(side_subvalues_list),
+ for j:1 thru length(zero_subintervals) do
+ (tmp5:delete(side_subvalues_list[zero_subintervals[j]],tmp5)),
+ non_zero_side_subvalues_list:tmp5,
+
+ /* if length(non_zero_side_subintervals_list)=0 then we are done just with token2 */
+
+ if is(equal(length(non_zero_side_subintervals_list),0)) then return(token2),
+
+ /* otherwise, we must also take into account token3 */
+ LL:length(non_zero_side_subintervals_list)/2,
+ for k:1 thru LL do count2[k]:charfun(is(equalp(non_zero_side_subvalues_list[k],ratsubst(-x,x,non_zero_side_subvalues_list[2*LL+1-k])))),
+ if is(equal(sum(count2[j],j,1,LL),LL)) then token3:0,
+ for k:1 thru LL do count3[k]:charfun(is(equalp(non_zero_side_subvalues_list[k],-ratsubst(-x,x,non_zero_side_subvalues_list[2*LL+1-k])))),
+ if is(equal(sum(count3[j],j,1,LL),LL)) then token3:1,
+ if is(not(equal(token3,0) or equal(token3,1))) then return(-1) elseif
+ is(token2#token3) then return(-1)
+ elseif is(equal(token2,0)) then return(0)
+ elseif is(equal(token2,1)) then return(1)
+
+ ),
+
+ /* we continue here with unbounded intervals */
+
+ expr1:subvalues[tmp2[1]],
+ expr2:ratsubst(-x,x,subvalues[tmp2[2]]),
+ if is(equalp(expr1,expr2)) then token1:0
+ elseif is(equalp(expr1,-expr2)) then token1:1
+ else return(-1),
+
+ /* the problem when defining token1 is that if the asymptotic value is 0, then it is always token1:0 */
+ /* although the function could be odd. To cope with this, we define token0 below */
+
+ tmp3:makelist(subintervals[i],i,sublist_indices(tmp1,lambda([x],x=0))),
+ if is(equal(length(tmp3),0)) then return(token1),
+ subvalues_list:makelist(subvalues[i],i,sublist_indices(tmp1,lambda([x],x=0))),
+ token0:if member(0,makelist(subvalues[i],i,tmp2)) then 0 else 7,
+ for j:1 thru length(tmp3) do
+ (
+ M[j]:makelist(part(tmp3[j],k),k,1,length(tmp3[j])),
+ for m:1 thru 2 do N[j,m]:makelist(part(M[j],m,n),n,1,length(M[j])),
+ P[j]:append(N[j,1],N[j,2]),Q[j]:sort(delete(x,P[j]),"<")
+ ),
+ subintervals_list:makelist(Q[j],j,1,length(tmp3)),
+ for i:1 thru length(subintervals_list) do (if is(lmin(subintervals_list[i])*lmax(subintervals_list[i])<0) then icentral:i),
+ if is(icentral>0) then middle:subvalues_list[icentral],
+ if is(icentral>0) then
+ (if evenfunp(middle,x) then token2:0 elseif oddfunp(middle,x) then token2:1 else return(-1))
+ else
+
+ /* now we analyze what happens if icentral=0, so there are only non-bounded intervals and side intervals */
+
+ (
+ /* as before, there must be an even number of symmetric intervals, otherwise -1 */
+ side_subintervals_list:copylist(subintervals_list),
+ if not(evenp(length(side_subintervals_list))) then return(-1),
+ L:length(side_subintervals_list)/2,
+ for k:1 thru L do count1[k]:charfun(is(equal(side_subintervals_list[k],reverse(map("-",side_subintervals_list[2*L+1-k]))))),
+ if sum(count1[j],j,1,L)#L then return(-1),
+ side_subvalues_list:copylist(subvalues_list),
+ for k:1 thru L do count2[k]:charfun(is(equalp(side_subvalues_list[k],ratsubst(-x,x,side_subvalues_list[2*L+1-k])))),
+ if is(equal(sum(count2[j],j,1,L),L)) then token3:0,
+ for k:1 thru L do count3[k]:charfun(is(equalp(side_subvalues_list[k],-ratsubst(-x,x,side_subvalues_list[2*L+1-k])))),
+ if is(equal(sum(count3[j],j,1,L),L)) then token3:1,
+ if is(not(equal(token3,0)) and not(equal(token3,1))) then return(-1) elseif
+ is(not(equal(token1,token3)) and not(equal(token0,0))) then return(-1)
+ elseif is( equal(token0,0) and equal(token3,0)) then return(0)
+ elseif is( equal(token0,0) and equal(token3,1)) then return(1)
+ elseif is(not(equal(token0,0)) and equal(token1,token3)) then return(token3)
+ elseif is(not(equal(token0,0)) and not(equal(token1,token3))) then return(-1)
+ ),
+
+ /* and what happens if icentral#0 */
+
+ side_subintervals_list:delete(subintervals_list[icentral],subintervals_list),
+ if is(equal(length(side_subintervals_list),0)) then (if is(equal(token1,token2)) then return(token1*token2) else return(-1)),
+ if not(evenp(length(side_subintervals_list))) then return(-1),
+ L:length(side_subintervals_list)/2,
+ for k:1 thru L do count1[k]:charfun(is(equal(side_subintervals_list[k],reverse(map("-",side_subintervals_list[2*L+1-k]))))),
+ if sum(count1[j],j,1,L)#L then return(-1),
+ side_subvalues_list:delete(subvalues_list[icentral],subvalues_list),
+
+ /* remove those subintervals in which funvar vanishes */
+
+ zero_subintervals:sublist_indices(side_subvalues_list,lambda([x],x=0)),
+ tmp4:copylist(side_subintervals_list),
+ for j:1 thru length(zero_subintervals) do
+ (tmp4:delete(side_subintervals_list[zero_subintervals[j]],tmp4)),
+ non_zero_side_subintervals_list:tmp4,
+ tmp5:copylist(side_subvalues_list),
+ for j:1 thru length(zero_subintervals) do
+ (tmp5:delete(side_subvalues_list[zero_subintervals[j]],tmp5)),
+ non_zero_side_subvalues_list:tmp5,
+
+ /* if length(non_zero_side_subintervals_list)=0 then we are done just with token1 and token2 */
+
+ if is(equal(length(non_zero_side_subintervals_list),0)) then (if is(equal(token1,token2)) then return(token1*token2) else return(-1)),
+
+ /* otherwise, we must also take into account token3 */
+ LL:length(non_zero_side_subintervals_list)/2,
+ for k:1 thru LL do count2[k]:charfun(is(equalp(non_zero_side_subvalues_list[k],ratsubst(-x,x,non_zero_side_subvalues_list[2*LL+1-k])))),
+ if is(equal(sum(count2[j],j,1,LL),LL)) then token3:0,
+ for k:1 thru LL do count3[k]:charfun(is(equalp(non_zero_side_subvalues_list[k],-ratsubst(-x,x,non_zero_side_subvalues_list[2*LL+1-k])))),
+ if is(equal(sum(count3[j],j,1,LL),LL)) then token3:1,
+ if is(not(equal(token0,0)) and is(not(equal(token1,token2)) or not(equal(token1,token3)) or not(equal(token2,token3))) ) then return(-1)
+ elseif is(not(equal(token0,0)) and is(equal(token1,0)) ) then return(0)
+ elseif is(not(equal(token0,0)) and is(equal(token1,1)) ) then return(1)
+ elseif is( equal(token0,0) and equal(token3,0) and equal(token2,0)) then return(0)
+ elseif is( equal(token0,0) and equal(token3,1) and equal(token2,1)) then return(1)
+ else return(-1)
+
+)$
+
+fourier_sec_coeff(fuvar,x):=
+block([pp,LL,lm,a0,coeff],
+ local(a,b,n,MM,NN,PP,QQ),
+ declare(n,integer),
+ pp:((length(fuvar)/2)-1),
+ LL:makelist(part(fuvar,i),i,makelist(2*s-1,s,1,pp)),
+ for j:1 thru length(LL) step 1 do
+ (
+ MM[j]:makelist(part(LL[j],r),r,1,length(LL[j])),
+ for r:1 thru 2 do NN(j,r):=makelist(part(MM[j],r,k),k,1,length(MM[j])),
+ PP[j]:append(NN(j,1),NN(j,2)),QQ[j]:delete(x,PP[j])
+ ),
+ for i:1 thru pp step 1 do partsums[i]:sort(QQ[i],"<"),
+ for i:1 thru pp step 1 do partfunc[i]:part(fuvar,2*i),
+ lm:lmax(unique(flatten(makelist(QQ[q],q,1,length(LL))))),
+ a0:(1/(2*lm))*sum(integrate(partfunc[i],x,partsums[i][1],partsums[i][2]),i,1,pp),
+ if is(equal(paritycheck(fuvar,x),1)) then a(n):=0 else a(n):=(1/lm)*sum(adefint(partfunc[i]*cos(%pi*n*x/lm),x,partsums[i][1],partsums[i][2]),i,1,pp),
+ if is(equal(paritycheck(fuvar,x),0)) then b(n):=0 else b(n):=(1/lm)*sum(adefint(partfunc[i]*sin(%pi*n*x/lm),x,partsums[i][1],partsums[i][2]),i,1,pp),
+ coeff:[a0,simplify_sum(a(n)),simplify_sum(b(n))],
+ factor(ratsimp(coeff))
+)$
+
+fourier_sec_series(fuvar,x,u):=
+block([pp,LL,lm,a0,coeff],
+ local(a,b,n,MM,NN,PP,QQ),
+ declare(n,integer),
+ pp:((length(fuvar)/2)-1),
+ LL:makelist(part(fuvar,i),i,makelist(2*s-1,s,1,pp)),
+ for j:1 thru length(LL) step 1 do
+ (
+ MM[j]:makelist(part(LL[j],r),r,1,length(LL[j])),
+ for r:1 thru 2 do NN(j,r):=makelist(part(MM[j],r,k),k,1,length(MM[j])),
+ PP[j]:append(NN(j,1),NN(j,2)),QQ[j]:delete(x,PP[j])
+ ),
+ for i:1 thru pp step 1 do partsums[i]:sort(QQ[i],"<"),
+ for i:1 thru pp step 1 do partfunc[i]:part(fuvar,2*i),
+ lm:lmax(unique(flatten(makelist(QQ[q],q,1,length(LL))))),
+ a0:(1/(2*lm))*sum(integrate(partfunc[i],x,partsums[i][1],partsums[i][2]),i,1,pp),
+ if is(equal(paritycheck(fuvar,x),1)) then a(n):=0 else a(n):=(1/lm)*sum(adefint(partfunc[i]*cos(%pi*n*x/lm),x,partsums[i][1],partsums[i][2]),i,1,pp),
+ if is(equal(paritycheck(fuvar,x),0)) then b(n):=0 else b(n):=(1/lm)*sum(adefint(partfunc[i]*sin(%pi*n*x/lm),x,partsums[i][1],partsums[i][2]),i,1,pp),
+ a0+intosum(sum(a(n)*cos(%pi*n*x/lm),n,1,u))+intosum(sum(b(n)*sin(%pi*n*x/lm),n,1,u)))$ \ No newline at end of file
diff --git a/texmf-dist/scripts/ketpic/ketlib/maximaL/matoperation.max b/texmf-dist/scripts/ketpic/ketlib/maximaL/matoperation.max
new file mode 100644
index 00000000..512a9ffb
--- /dev/null
+++ b/texmf-dist/scripts/ketpic/ketlib/maximaL/matoperation.max
@@ -0,0 +1,47 @@
+/* Maxima tips http://www.cymric.jp/maxima/tips.html */
+
+rowmult(X, n, c) := block([M],
+ M:matrix(),
+ for i:1 thru length(X) do
+ if i=n then M:addrow(M, c*row(X, i))
+ else M:addrow(M, row(X, i)),
+ M
+);
+rowadd(X, n, m, c) := block([M],
+ M:matrix(),
+ for i:1 thru length(X) do
+ if i=n then M:addrow(M, c * row(X, m) + row(X, i))
+ else M:addrow(M, row(X, i)),
+ M
+);
+rowch(X, m, n) := block([M],
+ M:matrix(),
+ for i:1 thru length(X) do
+ if i=m then M:addrow(M, row(X, n))
+ else if i=n then M:addrow(M, row(X, m))
+ else M:addrow(M, row(X, i)),
+ M
+);
+
+colmult(X, n, c) := block([M],
+ M:matrix(),
+ for i:1 thru length(transpose(X)) do
+ if i=n then M:addcol(M, c*col(X, i))
+ else M:addcol(M, col(X, i)),
+ M
+);
+coladd(X, n, m, c) := block([M],
+ M:matrix(),
+ for i:1 thru length(transpose(X)) do
+ if i=n then M:addcol(M, c * col(X, m) + col(X, i))
+ else M:addcol(M, col(X, i)),
+ M
+);
+colch(X, m, n) := block([M],
+ M:matrix(),
+ for i:1 thru length(transpose(X)) do
+ if i=m then M:addcol(M, col(X, n))
+ else if i=n then M:addcol(M, col(X, m))
+ else M:addcol(M, col(X, i)),
+ M
+);
diff --git a/texmf-dist/scripts/ketpic/ketlib/maximaL/poincare.mac b/texmf-dist/scripts/ketpic/ketlib/maximaL/poincare.mac
new file mode 100644
index 00000000..aefab3a2
--- /dev/null
+++ b/texmf-dist/scripts/ketpic/ketlib/maximaL/poincare.mac
@@ -0,0 +1,71 @@
+/* Poincare.mac is a set of routines for the computation
+of Poincare surfaces of sections of Hamiltonian systems */
+
+/* Author: Jose A Vallejo
+Universidad Autonoma de San Luis Potosi (Mexico)
+josanv@gmail.com */
+
+/* To avoid messages when loading draw, execute
+with_stdout("/dev/null",load(draw))$ (in Linux)
+with_stdout("NUL",load(draw))$ (in Windows) */
+
+/* Poincare.mac requires a compiled version of
+rkfun.lisp, a LISP implementation of the RK4
+algorithm due to Richard Fateman. The original
+code can be downloaded from
+https://people.eecs.berkeley.edu/~fateman/lisp/rkfun.lisp */
+
+/*load("rkfun.fasl");*/
+/*load("rkfun.lisp");*/
+
+
+hameqs(H,name):=block(
+[vv,t,tvv,n,Q,P,eqq,eqp,eqs],
+vv:args(lhs(apply(fundef,[H]))),
+tvv:cons(t,vv),
+n:length(vv)/2,
+Q:makelist(vv[2*j-1],j,1,n),
+P:makelist(vv[2*k],k,1,n),
+eqq:makelist(float(diff(apply(H,vv),P[j])),j,1,n),
+eqp:makelist(float(-diff(apply(H,vv),Q[j])),j,1,n),
+eqs:join(eqq,eqp),
+for j:1 thru 2*n do define(funmake(concat(name,j),tvv),block([],mode_identity(float,vv),eqs[j])),
+apply(compile,makelist(concat(name,j),j,1,2*n)),
+[makelist(apply(concat(name,j),tvv),j,1,2*n),vv,makelist(concat(name,j),j,1,2*n)]
+)$
+
+lextract(ll,n):=block([l,a,b],
+l:length(ll),
+a:rest(ll,n),
+b:rest(ll,n-l-1),
+append(b,a)
+)$
+
+poincare3d(H,name,inicond,timestep,coord):=block(
+[heq,vars,tvars,hfuns,c,sol],
+heq:hameqs(H,name),
+vars:heq[2],
+tvars:cons(t,vars),
+hfuns:heq[3],
+c:first(sublist_indices(vars,lambda([x],x=coord))),
+sol:rkfun(hfuns,tvars,float(inicond),timestep),
+map(lambda([x],lextract(x,c)),map(lambda([x],rest(x)),sol))
+)$
+
+poincare2d(H,name,inicond,timestep,scene):=block(
+[heq,vars,tvars,hfuns,solu,c,soluc,sola,solb,subind,sol,e,f],
+heq:hameqs(H,name),
+vars:heq[2],
+tvars:cons(t,vars),
+hfuns:heq[3],
+solu:rkfun(hfuns,tvars,float(inicond),timestep),
+c:first(sublist_indices(vars,lambda([x],x=first(scene)))),
+soluc:map(lambda([x],x[c+1]),solu),
+sola:rest(soluc,-1)-second(scene),
+solb:rest(soluc)-second(scene),
+subind:sublist_indices(sola*solb,lambda([x],is(x < 0))),
+sol:makelist(solu[k],k,subind),
+e:first(sublist_indices(vars,lambda([x],x=third(scene))))+1,
+f:first(sublist_indices(vars,lambda([x],x=fourth(scene))))+1,
+makelist([j[e],j[f]],j,sol)
+)$ \ No newline at end of file
diff --git a/texmf-dist/scripts/ketpic/ketlib/maximaL/rkfun.lisp b/texmf-dist/scripts/ketpic/ketlib/maximaL/rkfun.lisp
new file mode 100644
index 00000000..70aafc7f
--- /dev/null
+++ b/texmf-dist/scripts/ketpic/ketlib/maximaL/rkfun.lisp
@@ -0,0 +1,99 @@
+;; variation on runge-kutta (rk) routine in complex_dynamics.lisp by villate@fe.up.pt
+;; this one with changed by fateman@berkeley.edu
+
+#| rkfun may be much faster than rk because it calls compiled programs.
+ rkfun (this file) should also be compiled to speed it up.
+ how to use:
+
+
+ first, read and understand the documentation for rk().
+ next, define the ODE right hand sides as functions as shown below
+
+fun1(y,x,x0):= x0$
+fun2(y,x,x0):= block([], mode_declare([y, x0,x],float),(1-x^2)*x0-x);
+compile(fun1)$
+compile(fun2)$
+showtime:all$
+
+last(rkfun([fun1,fun2],[x,xdot],[0.0,0.6],[t,0,50,0.01]));
+
+/* which is, in one test, 18 times faster than */
+
+last(rk([xdot,(1-x^2)*xdot -x],[x,xdot],[0.0,0.6],[t,0,50,0.01]));
+
+|#
+
+(defun $rkfun (funs vars initial domain ;; taken from complex_dynamics.lisp
+ &aux d u funlist k1 k2 k3 k4 r1 r2 r3 traj r
+ (it (mapcar #'coerce-float (cddr domain))))
+ (unless ($listp funs) (setq funs `((mlist) ,funs)))
+ (unless ($listp initial) (setq initial `((mlist) ,initial)))
+ (unless ($listp vars) (setq vars `((mlist) ,vars)))
+ (dolist (var (cdr vars))
+ (unless (symbolp var)
+ (merror (intl:gettext "rk: variable name expected; found: ~M") var)))
+ (unless (symbolp (cadr domain))
+ (merror (intl:gettext "rk: variable name expected; found: ~M")
+ (cadr domain)))
+ (setq vars (append '((mlist)) (list (cadr domain)) (cdr vars)))
+ (setq r (append `(,(car it)) (mapcar #'coerce-float (cdr initial))))
+ (setq funlist (cdr funs))
+
+ (setq d (/ (- (cadr it) (car it)) (caddr it)))
+ (setq traj (list (cons '(mlist) r)))
+ (do ((m 1 (1+ m))) ((> m d))
+ (progn ;;; was ignore-errors
+ (setq k1 (mapcar #'(lambda (x) (apply x r)) funlist))
+ (setq r1 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* 0.5d0 (caddr it) x)) k1)))
+ (push (+ (car r) (/ (caddr it) 2)) r1)
+ (setq k2 (mapcar #'(lambda (x) (apply x r1)) funlist))
+ (setq r2 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* 0.5d0(caddr it) x)) k2)))
+ (push (+ (car r) (/ (caddr it) 2)) r2)
+ (setq k3 (mapcar #'(lambda (x) (apply x r2)) funlist))
+ (setq r3 (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (caddr it) x)) k3)))
+ (push (+ (car r) (caddr it)) r3)
+ (setq k4 (mapcar #'(lambda (x) (apply x r3)) funlist))
+ (setq u (map 'list #'+
+ (mapcar #'(lambda (x) (* #.(/ 1.0 6.0d0) x)) k1)
+ (mapcar #'(lambda (x) (* #.(/ 1.0 3.0d0) x)) k2)
+ (mapcar #'(lambda (x) (* #.(/ 1.0 3.0d0) x)) k3)
+ (mapcar #'(lambda (x) (* #.(/ 1.0 6.0d0) x)) k4)))
+ (setq r
+ (append
+ `(,(+ (car it) (* m (caddr it))))
+ (map 'list #'+ (cdr r) (mapcar #'(lambda (x) (* (caddr it) x)) u))))
+ (push (cons '(mlist) r) traj)))
+ (when (< (car r) (cadr it))
+ (let ((s (- (cadr it) (car r))) )
+ (declare(double-float s)(special s)(optimize (speed 3)(safety 0)))
+ (progn ;; was ignore-errors
+ (setq k1 (mapcar #'(lambda (x) (mapply x r nil)) funlist))
+ (setq r1 (map 'list #'+ (cdr r) (mapcar #'(lambda (x)
+ (declare(double-float x))
+ (* s 0.5d0 x)) k1)))
+ (push (+ (car r) (* 0.5d0 s)) r1)
+ (setq k2 (mapcar #'(lambda (x) (mapply x r1 nil)) funlist))
+ (setq r2 (map 'list #'+ (cdr r) (mapcar #'(lambda (x)(declare(double-float x))
+ (* s 0.5d0 x)) k2)))
+ (push (+ (car r) (* 0.5d0 s)) r2)
+ (setq k3 (mapcar #'(lambda (x) (mapply x r2 nil)) funlist))
+ (setq r3 (map 'list #'+ (cdr r) (mapcar #'(lambda (x)(declare(double-float x))
+ (* s x)) k3)))
+ (push (+ (car r) s) r3)
+ (setq k4 (mapcar #'(lambda (x) (mapply x r3 nil)) funlist))
+ (setq u (map 'list #'+
+ (mapcar #'(lambda (x)(declare(double-float x))
+ (* #.(/ 1.0d0 6.0d0) x)) k1)
+ (mapcar #'(lambda (x)(declare(double-float x))
+ (* #.(/ 1.0d0 3.0d0) x)) k2)
+ (mapcar #'(lambda (x) (declare(double-float x))
+ (* #.(/ 1.0d0 3.0d0) x)) k3)
+ (mapcar #'(lambda (x)(declare(double-float x))
+ (* #.(/ 1.0d0 6.0d0) x)) k4)))
+ (setq r
+ (append
+ `(,(cadr it))
+ (map 'list #'+ (cdr r) (mapcar #'(lambda (x)(declare(double-float x))
+ (* s x)) u))))
+ (push (cons '(mlist) r) traj))))
+ (cons '(mlist) (nreverse traj))) \ No newline at end of file